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Abstract 

Sesquinaries,  Magnetics  and  Atmospheres: 

Studies  of  the  Terrestrial  Moons  and  Exoplanets 

by 

Michael  Nayak 

The  surface  brightness  of  Deimos,  groove  patterns  on  Phobos,  crustal 
magnetic  anomalies  on  the  Moon  and  the  composition  of  exoplanet  atmospheres 
represent  some  of  the  most  interesting  and  puzzling  questions  in  planetary  science. 
Why  is  Deimos  significantly  brighter  and  smoother  than  its  partner  moon  Phobos? 
What  is  the  origin  of  the  crater  chain  “grooves”  on  Phobos?  Are  the  magnetic 
anomalies  in  the  lunar  South  Pole-Aitken  basin  a  remnant  of  the  basin’s  formation,  or 
do  they  owe  their  existence  to  a  primordial  period  of  lunar  dynamo  activity?  And 
finally,  as  visible  wavelength  telescopes  are  designed  and  tested  for  space-based 
exoplanet  detections,  can  we  use  observed  albedo  spectra  to  determine  radius,  gravity, 
cloud  pressure  heights  and  atmospheric  compositions  for  these  planets?  I  use 
dynamical  modeling,  magnetic  inversions  and  Markov  Chain  Monte  Carlo  retrievals 
to  address  these  questions.  Major  findings  include  1)  the  likelihood  of  isotropic 
redistribution  of  reaccreted  ejected  material  on  Deimos,  2)  the  creation  of 
hemispherical  catenae  from  the  creation  of  primary  craters  on  Phobos,  which  match 
the  locations  and  geomorphology  of  several  existing  grooves  well,  3)  the  first 
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directional  magnetic  survey  of  South  Pole-Aitken  basin  anomalies,  and  a  larger  than 
expected  diversity  in  recovered  paleopole  directions,  and  4)  the  critical  importance  of 
considering  the  effects  of  planet  phase  in  exoplanet  atmosphere  retrievals;  changing 
planet  phase,  when  combined  with  low  signal-to-noise  observations,  can  cause 
several  orders  of  magnitude  of  uncertainty  in  atmospheric  methane  composition  and 
cloud  pressure  height,  among  others. 
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Chapter  1 
Introduction 


The  exponential  increase  in  computational  capabilities,  in  the  realm  of  both 
individual  computing  and  supercomputing,  has  enabled  a  large  jump  in  simulation 
capability  for  planetary  scientists.  This  is  especially  evident  in  the  three  subfields 
explored  in  this  dissertation  -  ejecta  dynamics,  magnetic  inversions  and  broad-scale 
retrievals  of  gas  giant  atmospheric  properties  -  allowing  this  work  a  complexity  and 
level  of  insight  that  might  not  have  been  possible  even  ten  years  ago.  This  work  may 
be  divided  into  three  general  partitions.  In  order,  I  inquire  after  geophysical  puzzles 
related  to  the  moons  of  Mars  (sesquinaries.  Chapters  2  and  3),  magnetic  anomalies  on 
Earth’s  moon  (magnetics.  Chapter  4)  and  finally,  atmospheric  retrievals  for  inferring 
properties  of  gas  giant  planets  beyond  our  solar  system  (exoplanets.  Chapter  5). 

The  first  puzzle  explored  here  concerns  the  physical  appearance  and 
geological  features  noted  on  both  of  Mars’  tiny  moons,  Phobos  and  Deimos.  Deimos 
is  smoother  and  brighter  than  its  inner  counterpart,  while  Phobos  bears  several 
mysterious  “groove” -like  marks  on  its  surface.  The  dichotomy  in  appearance  and 
groove  structure,  while  first  noted  in  the  Viking  lander  era  almost  fifty  years  ago,  has 
to  date  not  been  fully  resolved.  By  using  multiple  simulations  with  thousands  of  test 
particles,  I  make  the  case  that  “sesquinaries”  may  have  a  significant  role  to  play.  In 
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Chapter  2,  I  consider  Mars’  outer  moon  Deimos.  By  modeling  sesquinary  impact 
cratering  from  the  impact  that  formed  its  largest  crater,  Voltaire,  it  is  found  that 
sesquinary  ejecta  globally  resurface  Deimos  near-isotropically,  erasing  the  previous 
geological  record.  I  conclude  that  dating  the  surface  of  Deimos  is  likely  more 
challenging  than  previously  suspected.  In  Chapter  3,  I  consider  Mars’  inner  moon 
Phobos.  Simulations  of  just-escaping  sesquinary  ejecta  show  persistent  reaccretions 
in  low- velocity  chain-like  clusters  similar  to  catenae  on  Ganymede  and  Callisto.  The 
morphological  similarity  to  linear  pitted  chains  on  Phobos  suggests  a  link  to  the  long- 
debated  mysterious  grooves  on  Phobos;  I  conclude  that  these  catenae  present  the 
missing  piece  to  families  of  grooves  that  do  not  fit  well  to  a  tidal  model  for  groove 
origin  as  Phobos  spirals  toward  Mars. 

The  second  puzzle  explored  in  this  work  pivots  to  magnetics,  but  still  with  a 
sesquinary  flavor.  Chapter  4  considers  Earth’s  Moon.  Despite  having  no  global 
magnetic  field,  the  Moon  still  exhibits  anomalous  localized  crustal  magnetic  fields. 
The  origin  of  these  fields  is  still  mysterious,  particularly  in  the  large  lunar  South  Pole- 
Aitken  basin.  Inverse  regression  techniques  applied  to  magnetic  anomalies  in  this 
region  show  diverse  paleopoles  that  imply  geophysically  improbable  amounts  of  true 
polar  wander.  A  number  of  possible  formation  hypotheses  are  explored,  including 
long  and  short  timescale  true  polar  wander  and  subsurface  dikes  magnetized  during 
an  ancient  period  of  lunar  dynamo  activity.  Continuing  the  thread  from  Chapters  2 
and  3,  it  is  found  that  secondary  ejecta  from  sesquinary  impactors,  reaccreting  as  the 
Moon  experienced  large  changes  in  its  moments  of  inertia  in  the  aftermath  of  the 
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South  Pole  impact,  might  also  be  another  plausible  explaination  for  the  large 
paleopole  variations  observed. 

For  my  fourth  and  final  study  (Chapter  5),  I  pivot  to  the  detection  of 
extrasolar  planets  (exoplanets).  The  atmospheric  and  physical  properties  of  these 
gas  giants  are  explored;  I  use  parallelized  Markov  Chain  Monte  Carlo  (MCMC) 
techniques  to  perform  atmospheric  retrievals  on  future  targets  of  optical  wavelength 
imaging,  and  retrieve  properties  such  as  methane  abundance,  planet  radius,  planet 
gravity  and  cloud  properties.  The  inverse  recovery  of  spectral  signatures  of  exoplanet 
coronagraph  targets,  in  the  visible  wavelength  regimes  explored  here,  bear  particular 
application  to  the  planning  phase  of  an  upcoming  NASA  direct-imaging  exoplanet 
mission  (WFIRST). 

In  summary,  this  dissertation  aims  to  i)  expand  studies  of  sesquinary  impacts 
and  help  quantify  their  geophysical  importance;  2)  present  implications  inferred  from 
diverse  paleopoles  from  lunar  magnetic  anomalies,  and  3)  implement  retrieval 
techniques  to  study  the  atmospheres  of  gas  giant  planets  in  solar  systems  beyond  our 
own. 
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Chapter  2 

Sesquinaries:  A  Case  Study  of  Deimos 


This  chapter  is  a  modified  reprint  of  M.  Nayak,  F.  Nimmo  and  B.  Udrea  (2016), 
Effects  of  Mass  Transfer  between  Martian  Satellites  on  Surface  Geology,  Icarus  267, 
pp.  220-231,  DOI:  10.1016/j.icarus.2015. 12.026. 

2.1  Abstract 

Impacts  on  planetary  bodies  can  lead  to  both  prompt  secondary  craters  and 
projectiles  that  reimpact  the  target  body  or  nearby  companions  after  an  extended 
period,  producing  so-called  “sesquinary”  craters.  This  chapter  examines  sesquinary 
cratering  on  the  moons  of  Mars.  By  modeling  the  impact  that  formed  Voltaire,  the 
largest  crater  on  the  surface  of  Deimos,  the  orbital  evolution  of  resulting  high- 
velocity  ejecta  across  500  years  is  explored  using  four-body  physics  and  particle 
tracking. 

The  bulk  of  mass  transfer  to  Phobos  occurs  in  the  first  10^  years  after  impact, 
while  reaccretion  of  ejecta  to  Deimos  is  predicted  to  continue  out  to  a  10"^  year 
timescale,  in  agreement  with  [Soter,  1971].  Relative  orbital  geometry  between 
Phobos  and  Deimos  plays  a  significant  role;  depending  on  the  relative  true  longitude, 
mass  transfer  between  the  moons  can  change  by  a  factor  of  five.  Of  the  ejecta  with  a 
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velocity  range  capable  of  reaching  Phobos,  25  -  42%  by  mass  reaccretes  to  Deimos 
and  12  -  21%  impacts  Phobos.  Ejecta  mass  transferred  to  Mars  is  <10%. 

It  is  found  here  that  the  characteristic  impact  velocity  of  sesquinaries  on 
Deimos  is  an  order  of  magnitude  smaller  than  those  of  background  (heliocentric) 
hypervelocity  impactors  and  will  likely  result  in  different  crater  morphologies.  The 
time-averaged  flux  of  Deimos  material  to  Phobos  can  be  as  high  as  11%  of  the 
background  (heliocentric)  direct-to-Phobos  impactor  flux.  This  relatively  minor 
contribution  suggests  that  spectrally  red  terrain  on  Phobos  [Murchie  and  Erard,  1996] 
is  not  caused  by  Deimos  material.  However  the  high-velocity  ejecta  mass  reaccreted 
to  Deimos  from  a  Voltaire-sized  impact  is  comparable  to  the  expected  background 
mass  accumulated  on  Deimos  between  Voltaire-size  events.  Considering  that  the 
high-velocity  ejecta  contains  only  0.5%  of  the  total  mass  sent  into  orbit,  sesquinary 
ejecta  from  a  Voltaire-sized  impact  could  feasibly  resurface  large  parts  of  the  moon, 
erasing  the  previous  geological  record.  Dating  the  surface  of  Deimos  may  be  more 
challenging  than  previously  suspected. 

2.2  Introduction 

Several  features  about  the  surface  geology  on  the  moons  of  Mars  remain 
poorly  understood.  The  grooves  on  Phobos,  which  do  not  exist  on  Deimos,  have 
received  the  most  attention  [Thomas,  1979;  Weidenschilling,  1979;  Horstman  and 
Melosh,  1989],  and  theories  for  their  formation  continue  to  be  proposed  to  this  day 
[Murray  et  al,  2006;  Hamelin,  2011;  Basilevsky  et  al,  2014;  Asphaug  et  ah,  2015b; 
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Nayak  and  Asphaug,  2015,  2016;  Wilson  and  Head,  2015].  However  this  is  far  from 
the  only  mystery.  Though  both  moons  are  heavily  cratered,  with  saturated  surfaces 
and  fine-grained  regolith  from  impact  debris  accumulation  [Thomas,  1979;  Lunine  et 
al,  1982],  a  large  portion  of  ejecta  produced  on  Deimos  is  retained  in  the  form  of 
crater  fill  of  ~5  m  depth  [Thomas  and  Veverka,  1980b],  a  phenomenon  not  noted  on 
Phobos.  This  difference  is  still  unexplained  [Lee,  2009].  The  surface  of  Deimos  is 
also  significantly  smoother  and  brighter  than  Phobos,  likely  a  result  of  crater  fill 
[Veverka,  1978;  Thomas,  1993;  Thomas  et  al.,  1996]. 

Phobos  also  exhibits  two  distinct  spectral  units,  one  of  “redder”  origin  and  one 
of  “bluer”  origin,  possibly  stemming  from  a  compositional  difference  [Murchie  and 
Erard,  1996;  Lee,  2009].  The  bluer  unit  is  associated  with  the  Stickney  crater  and  an 
origin  from  depth.  The  redder  unit  associated  with  the  surface  and  small  craters;  it  is 
spectrally  similar  to  D-type  asteroids,  but  also  to  Deimos  [Murchie  and  Erard,  1993, 
1996].  It  has  been  proposed  that  the  red  unit  is  a  wide-spread  shallow  layer 
superimposed  on  a  blue  base  [Murchie  and  Erard,  1996],  for  which  there  are  four 
possible  causes  [Britt  and  Pieters,  1988;  Murchie  et  al,  1991]:  1)  accretion  of  D- 
asteroid  material  onto  blue  Phobos  material;  2)  optical  alteration  of  the  bluer  unit;  3) 
accretion  of  ejecta  from  Martian  basin  impacts  and  subsequent  space  weathering  ox  4) 
Phobos  is  an  inherently  heterogeneous  rubble  pile  and  the  red/blue  units  are  end- 
member  compositions.  One  aim  of  this  study  is  to  investigate  the  possibility  that  the 
red  veneer  on  top  of  the  base  blue  unit  may  be  ejecta  accreted  from  Deimos  rather 
than  Mars. 
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Previous  work  has  established  that  impact  ejecta  can  reimpact  the  target  body 
or  nearby  companions  after  an  extended  period,  creating  so-called  “sesquinary” 
impact  morphology.  Examples  of  sesquinary  studies  in  the  literature  include  Earth’s 
Moon  [Gladman  et  ah,  1995],  lo  [Alvarellos  et  ah,  2008],  Ganymede  [Alvarellos  et 
al,  2002],  Europa  [Zahnle  et  al,  2008]  and  Pluto  [Bierhaus  and  Danes,  2014],  Eor 
Mars,  previous  work  suggests  ejecta  released  at  slightly  greater  than  the  satellite’s 
escape  velocity  could  remain  in  the  system  and  subsequently  reimpact  at  low  relative 
velocities  [Safer,  1971,  1972],  Possible  evidence  for  this  was  noted  in  analysis  of 
Viking  images  [Veverka  and  Duxbury,  1977],  however  the  efficiency  of  this  process 
was  previously  unknown  [Thamas,  1979].  This  chapter  reports  on  the  distribution  of 
impact  velocities  and  geometries  from  inter-moon  mass  transfer  trajectories,  and 
present  conclusions  on  the  role  and  importance  of  sesquinary  mass  transfer  between 
the  Martian  moons. 

2.3  Methods 

Voltaire,  the  largest  confirmed  crater  on  Deimos,  has  a  diameter  of  3  km 
[Thamas,  1979;  Thamas  and  Veverka,  1980a].  By  modeling  the  orbital  evolution  of 
ejecta  from  the  Voltaire-forming  impact,  we  aim  to  characterize  an  end-member  case 
of  mass  transfer  from  Deimos  to  other  Martian  system  bodies. 

To  model  the  streamlines  ejected  by  the  Voltaire  impact,  we  use  a  simplified 
form  of  Maxwell’s  Z-model  [Maxwell  and  Seifert,  1974;  Maxwell,  1977;  Raddy, 
1977].  Eirst,  coordinate  transformations  necessary  to  use  surface-centered  Z-model 
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streamlines  in  a  Mars-centered  simulation  are  detailed.  Three  planet-based  frames  are 
used,  explained  below. 

2.3.1  Coordinate  Transformations 

2.3.1. 1  Planet  Centered  Inertial  (PCI)  frame 

The  origin  of  the  PCI  frame  is  the  center  of  the  body:  Mars,  Deimos  or 
Phobos.  The  positive  x-axis  points  toward  the  vernal  equinox,  the  positive  z-axis 
extends  through  the  North  Pole  of  the  planet  and  the  y-axis  completes  the  right  hand 
system.  In  this  definition,  the  North  Pole  is  that  pole  of  rotation  that  lies  on  the  north 
side  of  the  invariable  plane  of  the  solar  system  [Archinal  etal,  2010].  The  planetary 
system  model  described  in  Section  2.3.4  is  placed  with  reference  to  the  Mars 
Centered  Inertial  (MCI)  frame. 

2.3.1.2  Planet  Centered  Planet  Fixed  (PCPF)  frame 

Like  the  PCI  frame,  the  origin  of  the  PCPF  frame  is  also  the  center  of  the 
body,  and  shares  its  z-axis  definition.  However,  the  x-axis  extends  through  the 
intersection  between  the  planet’s  equator  and  its  prime  meridian,  with  the  y-axis 
completing  the  right  hand  system. 

Figure  2.1  illustrates  the  relationship  between  the  PCI  and  PCPF  frames.  PCI 
can  be  rotated  into  the  PCPF  frame  around  the  z-axis  with  the  rotation  matrix: 

r  cosw  sinw  O' 

PCPFj^pci  —  —sinw  cosw  0  (2.1) 

.0  0  1. 

where  w  is  the  angle  of  rotation. 
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2.3. 1.3  Topocentric  Horizon  Frame  (also  SEZ:  South-East-Zenith  frame) 

The  topocentric  horizon  frame  is  adapted  from  the  South-East-Zenith  (SEZ) 
frame  as  defined  by  [Vallado,  2013];  the  two  are  referred  to  interchangeably  here. 
The  definition  of  the  Topocentric  Horizontal  frame  assumes  a  sphere  centered  at  the 
center  of  mass  of  Deimos  and  tangent  to  the  origin  of  the  frame,  which  is  the  center  of 
the  impact  site  (Voltaire).  The  x-axis  is  aligned  with  the  meridian  that  passes  through 
the  center  of  Voltaire  and  points  south.  The  y-axis  is  defined  such  that  the  x-y  plane  is 
tangent  to  the  surface  of  Deimos  at  the  center  of  Voltaire  and  points  along  the  local 
latitude  circle.  Completing  the  right-handed  system,  the  z-axis  points  radially  outward 
from  the  impact  site  towards  the  “local”  zenith.  The  local  horizon  forms  the 
fundamental  plane  for  this  system,  i.e.,  the  plane  defined  by  the  south  and  east  axes.  It 
should  be  noted  that  there  is  a  subtle  difference  between  the  definition  of  the  impact 
site’s  latitude  by  geodetic  or  astronomical  standards  [Vallado,  2013];  these  become 
identical  by  imposing  the  assumption  of  a  perfectly  spherical  impacted  body 
(Deimos).  The  low  gravitational  acceleration  at  the  surface  of  Deimos  has  a 
negligible  effect  on  the  speed  of  the  ejecta,  so  the  assumption  of  a  uniform  spherical 
geometry  is  justified. 

Eigure  2.2  illustrates  the  relationship  between  SEZ  and  PCPE.  is  the 
elevation  angle  of  the  ejection  velocity  vector  from  the  horizontal,  defined  as 
0  <  £gj  <  90°.  Prom  the  Z-model  formulation  (see  Section  2.3.2),  Sgj  =  35.4°.  Pgj  is 
the  azimuth  of  the  ejection  velocity  vector  and  is  measured  from  the  North,  positive 
clockwise  as  viewed  from  above  the  site  such  that  0  <  <  360°,  and  is  sampled 
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at  1 1  positions  across  this  range  separated  by  30°.  The  streamlines  defined  in  Z-model 
frame  can  be  rotated  into  SEZ  using  the  relationships: 

Vsouth  =  Vr  COS(180°  -  pej) 

VEast  =  t^rSin  (180°  -  pej)  (  2.2  ) 

where  Vr  is  defined  by  Equation  2.7.  Subsequently  two  rotation  matrices,  the 
first  around  the  y-axis  and  the  second  around  the  z-axis,  are  required  to  rotate  the  SEZ 
frame  into  the  PCPE  frame  (specifically,  the  Voltaire  SEZ  frame  into  the  Deimos- 
Centered  Deimos-Eixed  frame).  The  rotation  matrices  are: 


cos  Aj 

—sin  Aj 

O' 

-  sin  (Pi 

0 

cos  (p{ 

DCDFj^SEZ^  SinA^ 

cos  Aj 

0 

0 

1 

0 

-  0 

0 

1- 

cos  (Pi 

0 

sin  (Pi. 

where  (pi  and  Aj  are  illustrated  in  Eigure  2.2.  Einally,  the  DCDE  velocity 
coordinates  are  converted  to  MCI.  By  manipulating  the  basic  kinematic  equation  for 
the  position  vector  of  an  ejected  particle  in  the  DCDE  frame,  in  terms  of  the 

known 

^MCl  ^DCDF  +  MCI^DCDF  ^  ^DCDF  (  2.4  ) 

Now  ,  where  v^^'is  the  velocity  vector  of  Deimos  in  the  MCI 

frame.  This  is  determined  from  Deimos  ephemerides.  jg  angular  rate  of 

the  DCDE  frame  with  respect  to  MCI;  this  is  the  rotation  rate  of  Deimos.  is  the 

position  of  the  impact  site  in  DCDE  coordinates. 
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Following  a  similar  process,  for  the  impact  of  ejecta  with  Phobos  the  process 
is  reversed  to  obtain  the  velocity  vector  at  impact  in  the  Phobos  Centered  Phobos- 
Fixed  (PhCPhF)  frame.  With  similar  notation  as  used  above  this  relation  is: 

.-IfPhCPhF  _  .-liMCI  MCI  .PhCPhF  w  ;^PhCPhF  I  1  K  \ 

^ej  ~  ^ej  ^Ph  CO  X  z.o  ) 

The  validity  of  these  equations  has  been  checked  with  MICE,  a  commercial 
level  interface  created  by  JPL/Caltech  to  SPICE  ephemeris  information  from  NASA’s 
Navigation  and  Ancillary  Information  Eacility  (naif.jpl.nasa.gov)  [Acton  et  al,  2002]. 

2.3.2  Impact  Model:  Generating  2-D  Streamlines 

To  model  the  streamlines  ejected  by  the  Voltaire  impact,  we  use  a  simplified 
form  of  Maxwell’s  Z-model  [Maxwell  and  Seifert,  1974;  Maxwell,  1977;  Roddy, 
1977].  Though  limited  by  its  neglect  of  interactions  across  streamlines,  the  Z-model 
reasonably  approximates  several  experimentally  observed  features  [Melosh,  1989; 
Richardson  et  al,  2007].  The  limitations  of  a  Z-model  implementation  are  discussed 
at  length  by  [Barnhart  and  Nimmo,  2011].  This  application  is  only  concerned  with 
ejecta  streamlines  that  escape  Deimos,  and  is  unaffected  by  the  details  of  cratering 
flow  beneath  the  ground  plane,  surface  material  mixing  during  ejection  or  direct 
retention  and  emplacement  of  deposits.  Therefore  it  provides  a  suitable  level  of 
insight  into  an  outbound  velocity  distribution;  approximations  made  by  the  Z-model 
are  unlikely  to  alter  our  qualitative  results. 

The  formulation  of  [Barnhart  and  Nimmo,  201 1]  is  adopted  here,  who  use  Z  = 

2.71  for  a  Mars  application.  When  tested  against  numerical  computations,  Z  =  2.7 

represents  surface  explosion  cratering  flow  well  [Melosh,  1989].  All  streamlines  are 

11 


ejected  at  a  constant  angle  of  35.4°  from  the  horizontal,  set  according  to  the  relation 
[Maxwell,  1977]: 

Eej  =  tan“^(Z  —  2)  (  2.6  ) 


Outbound  radial  (Vr)  and  vertical  (Vz)  ejection  velocities  vary  inversely  with 
distance  from  the  center  of  the  crater  r  [Maxwell,  1977]  as: 

Vr  =  a/r^  (  2.7  ) 

Vz  =  (Z-  2)Vr  (  2.8  ) 

where  go  denotes  the  acceleration  due  to  gravity  for  Deimos  (0.003  m/s^)  and: 


a  = 


Qo^t 


2Z+1 


4Z(Z-2) 


(2.9) 


Using  a  final  crater  radius  R/  =  1500  m  for  Voltaire  [Veverka,  1978],  the 
transient  crater  radius  is  calculated  as  R?  =  0.65  R/  [Barnhart  and  Nimmo,  201 1].  For 
the  analysis  presented  here  the  number  of  streamlines  (n)  has  been  chosen  to  yield  a 
suitably  dense  streamline  distribution  with  velocities  greater  than  the  Deimos  escape 
velocity.  Setting  Rmin=  0  varying  Rmin  <  r  <  R(,  n  =  600  streamlines  evenly 
spaced  in  radius  are  generated  within  the  Voltaire  crater.  Converting  streamlines  into 
axisymmetric  coordinates  [cf.  Barnhart  and  Nimmo,  2011,  Fig  1],  the  radial  and 
vertical  coordinates  are  extracted  as: 

r  =  Rj  sin  0(1  —  cos  0)^  (2.10) 

z  =  Rj  cos  0(1  —  cos  0)^  (2.11) 

where  0  is  the  angle  from  the  vertical  { 0  |  0  G  0  :  7r/2 }  and: 

^  «r±nm  2.12) 

‘  n 
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2.3.3  Creating  3-D  Velocity  Streamlines 


The  2-D  axisymmetrical  distribution  is  now  used  to  create  an  approximation 
to  a  3-D  excavation.  The  fate  of  the  ejecta  particle  (reaccretion  to  Deimos,  impact  to 
Phobos,  impact  to  Mars  or  escape)  can  vary  greatly  depending  on  the  azimuth  of  the 
streamline.  To  rotate  around  the  azimuthal  direction,  we  define  the  Topocentric 
Horizon  frame  (Section  2. 3. 1.3),  adapted  from  the  South-East-Zenith  (SEZ)  frame 
[Vallado,  2013].  The  azimuth  of  the  ejection  velocity  vector  P  is  measured  from  the 
north,  clockwise  as  viewed  from  above  the  impact  site.  A  30°  span  is  selected  as  a 
compromise  between  computational  efficiency  and  sampling  a  variety  of  azimuths 
across  the  possible  solution  space,  such  that  (^  \  P  E(0:  n/6:  2n)  for  a  total  of  11 
possible  azimuths.  This  yields  a  three-dimensional  outbound  velocity  distribution  tied 
to  Voltaire.  Eor  use  with  the  Mars  gravity  system  integrator,  these  coordinates  are 
then  rotated  into  the  Mars  Centered  Inertial  (MCI)  frame;  details  of  coordinate 
transformations  through  the  Deimos-Centered  Deimos-Eixed  (DCDE)  and  Deimos- 
Centered  Inertial  (DCI)  frames  were  presented  in  Section  2. 3. 1.3. 

Einally,  this  work  is  specifically  interested  in  those  streamlines  that  have 
sufficient  velocity  to  reach  the  orbit  of  Phobos.  Since  both  moons  lie  in  the  same 
orbital  plane  \Cazenave  et  al,  1980],  the  minimum  velocity  at  Deimos  to  reach 
Phobos  can  be  analytically  calculated  with  the  Hohmann  transfer  [Section  6.3,  Curtis, 
2013].  Particles  begin  to  cross  the  orbit  of  Phobos  at  velocities  above  500  m/s,  so  we 
set  the  lower  bound  on  velocities  of  interest  at  400  m/s.  Prom  Deimos,  the  minimum 
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velocity  to  escape  the  gravitational  well  of  Mars  is  analytically  approximated  as  [Eqn 
2.80,  Curtis,  2013]: 

^esc  "^i^/^Deimos  (  2.13  ) 

where  VDeimos  is  the  distance  from  Deimos  to  Mars  and  jU  is  the  product  of  the 
gravitational  constant  and  the  mass  of  Mars.  From  Equation  2.13,  Vgsc  =  1.91  km/s; 
we  set  the  upper  bound  on  velocities  of  interest  at  2  km/s.  Therefore,  velocity 
streamlines  in  the  range  {v  \  v  E  400  :  2000  m/s}  are  examined.  Nineteen  of  600 
streamlines  fall  within  this  range;  rotated  around  1 1  azimuthal  positions,  this  creates  a 
209-streamline  distribution.  While  this  work  focuses  on  ejecta  with  sufficient  velocity 
to  reach  Phobos  (-400  m/s),  note  that  the  majority  of  ejecta  launched  from  Deimos  at 
lesser  velocities  will  ultimately  re-impact  Deimos. 

2.3.4  Planetary  System  Model 

This  section  details  the  formulation  of  the  planetary  model.  Centered  at  the 
primary,  the  Mars  gravity  system  is  modeled  with  12x12  gravity  harmonics  from  the 
NASA  Planetary  Data  System  [pds-geosciences.wustl.edu]  [Murchie,  2010].  The 
effects  of  permanent  solid  tides  are  included,  truncated  to  the  size  of  the  gravity  field. 
The  present-day  orbit  of  Deimos  is  likely  similar  to  its  primordial  orbit  [Burns,  1978; 
Lambeck,  1979];  the  500-year  orbits  of  Phobos  and  Deimos  are  generated  analytically 
from  modern-day  mean  orbital  parameters  [ssd.jpl.nasa.gov.  Table  2.1].  A  subset  of 
streamlines  was  run  against  high-precision  orbits  generated  for  Phobos  and  Deimos 
[Genova  and  Folkner,  personal  communication,  2075];  results  were  not  found  to 
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differ  substantially  from  those  run  against  the  analytical  orbits.  In  the  interest  of 
computational  speed,  the  analytical  formulation  is  adopted  hereafter. 

Table  2.1.  Mean  Orbital  Parameters  and  Constants  for  Phobos  and  Deimos 


Deimos 

Phobos 

Semi-major  axis 

23485  km 

9389.8  km 

Eccentricity 

0.00115571 

0.0164255 

Inclination 

1.79  deg 

1.09  deg 

Right  Ascension  of  Ascending  Node 

148.0  deg 

319.9  deg 

Argument  of  Periapsis 

123.3  deg 

270.7  deg 

Mean  longitude^ 

109.7  deg 

190.6  deg 

Rate  of  mean  longitude 

0.00330049  deg/s 

0.0130317  deg/s 

Mean  radius 

6.2  km 

11.3  km 

Acceleration  due  to  gravity 

0.003  m/s^ 

0.0057  m/s^ 

Escape  velocity 

5.56  m/s 

11.39  m/s 

Hill  sphere  radius 

16.5  km 

21.5  km 

Due  to  its  proximity  to  Mars,  the  orbit  of  Deimos  is  primarily  influenced  by 
Mars’  oblateness;  the  third-body  effect  from  the  Sun  or  other  planetary  bodies  such  as 
Jupiter  is  negligible  [Burns,  1972].  Ejecta  released  from  the  orbits  of  Deimos  will 
follow  a  similar  pattern;  we  therefore  neglect  these  third-body  effects.  Similarly,  solar 
radiation  pressure  is  a  second-order  effect  when  compared  to  solar  gravity 
perturbations  [Klacka,  2002;  Farnocchia  et  al,  2014];  we  neglect  this  effect  as  well. 
However  for  complete  understanding  of  orbital  evolution  within  the  Martian  system 
we  include  third-body  perturbation  effects  from  Phobos  and  Deimos,  making  the 
physics  of  our  model  a  four-body  problem. 

It  is  assumed  that  any  particle  that  enters  the  Hill  sphere  (Table  2.1)  of  either 
moon  will  be  captured  by  it^.  Due  to  its  irregularly  triaxial  shape,  Deimos  has  an 


^  Mean  longitude  is  calculated  with  reference  to  a  mean  epoch  coordinate  system:  The  mean  equator-mean  equinox 
coordinate  system  is  evaluated  at  the  epoch  of  the  object.  The  starting  epoch  is  arbitrary  due  to  our  evaluation  of  the 
orbital  dynamics  at  multiple  relative  geometry  configurations  between  Mars,  Phobos  and  Deimos  that  encompass  all 
possible  geometries  between  the  bodies. 
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uneven  gravity  field  that  causes  the  escape  velocity  to  be  lower  at  the  sub-Mars  and 
anti-Mars  points  [Davis  et  al,  1981].  Ejecta  in  the  4-6  m/s  range  will  see  the  largest 
variation  in  range  [Thomas,  1993];  since  the  slowest  particle  we  consider  is  ejected  at 
>400  m/s  and  the  escape  velocity  varies  on  the  order  of  cm/s,  it  may  be  safely 
assumed  that  the  escape  velocity  at  Voltaire  equals  the  average  escape  speed  over 
Deimos. 

Finally  the  impact  of  relative  orbital  geometry  is  considered.  Though  an 
analytical  formulation  has  been  used  to  consider  similar  problems  in  the  past  [Soter, 
1971;  Dobrovolskis  and  Bums,  1980;  Thomas,  1998],  this  approach  may  be 
insufficient  for  a  full  understanding  of  ejecta  dynamics.  Phobos  is  closer  to  Mars  than 
any  other  planetary  satellite,  and  is  the  only  moon  with  an  orbital  period  less  than  the 
rotational  period  of  its  primary  body  [Burns,  1972].  The  flux  of  material  impacting 
Phobos  can  vary  drastically  between  inferior  and  superior  conjunctions  between 
Phobos  and  Deimos.  The  difference  in  the  true  longitude  between  Deimos  and 
Phobos  can  (and  does,  see  Figure  2.3)  change  the  outcome  of  a  Phobos  collision  to  a 
Mars  collision,  or  vice  versa.  Therefore,  though  the  orbits  of  Phobos  and  Deimos  are 
generated  analytically,  all  propagation  in  this  work  is  ephemeris-centered. 

Deimos  has  an  orbital  period  of  30.3  hours,  and  Phobos  7.5  hours.  To  evaluate 
the  fate  of  ejecta  across  the  range  of  possible  Mars-Deimos  orbital  geometries,  this 
orbit  is  discretized  into  28  geometry  configurations  (GCs),  evenly  spaced  in  one-hour 

n 

To  test  the  validity  of  using  the  Hill  sphere  as  an  impact  boundary,  we  selected  100  Phobos  impact  trajectories  at 
random  and  integrated  them  with  an  impact  boundary  of  13  km,  the  longest  semi-major  axis  of  Phobos  [Murchie  et  al, 
2003];  98/100  trajectories  were  still  found  to  impact.  Therefore,  our  qualitative  results  are  not  affected  by  the  use  of  the 
Hill  sphere  as  an  impact  boundary. 
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increments  from  the  Deimos  apoapsis.  In  this  time  period,  Phobos  completes  nearly 
four  orbits  of  Mars,  allowing  for  discretization  of  the  range  of  possible  Mars-Phobos- 
Deimos  orbital  geometries  as  well.  At  each  GC,  209  streamlines  are  released  and 
propagated,  for  a  total  of  5,852  streamlines,  thereby  ensuring  a  robust  capture  of  the 
impacting  process  despite  variations  in  orbital  positions  and  conjunction  geometries. 

Estimates  for  lifetime  of  ejecta  in  the  Martian  system  range  from  10  to  10 
years  [Soter,  1971;  Davis  et  al,  1981].  Each  streamline  is  integrated  in  the  Mars 
gravity  system  for  tmax  =  500  years,  stopping  sooner  only  in  the  event  of  a  planetary 
body  collision  or  departure  from  the  Mars  gravitational  sphere  of  influence.  This 
length  of  integration  balances  computational  feasibility  with  permitting  a  statistically 
significant  number  of  ejected  particles  to  impact  or  escape.  As  shall  be  shown,  the 
uncertainty  introduced  by  doing  so  does  not  affect  our  conclusions.  It  also  permits  the 
use  of  a  Runge-Kutta  integrator  without  excessive  approximations  to  the  perturbed 
Hamiltonian  [Leimkuhler  and  Reich,  2004].  A  seventh-order  Runge-Kutta  integrator 
with  eighth-order  error  control  is  used  for  all  orbit  propagations.  The  maximum 
permitted  relative  error  is  10  '°.  The  Tisserand  parameter  is  used  to  evaluate  the 
performance  of  the  integrator,  according  to  which: 

^  -I-  2^/ a(l  —  e^)cos  i  =  constant  (2.14) 

where  a,  e  and  i  are  the  semi-major  axis,  eccentricity  and  inclination  of  the 
orbit.  The  differences  in  the  Tisserand  parameter  for  an  individual  particle  are  no 
greater  than  10'^,  i.e.,  at  most,  a  0.001%  change  across  the  integration  period. 
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2.4  Results 

Figure  2.3  plots  the  fate  of  5,852  massless  ejecta  particles  in  the  Martian 
gravity  system  across  500  years  as  a  function  of  the  orbital  geometry  and  Voltaire 
ejection  velocity.  The  number  of  reaccretions  to  Deimos  and  particles  still  flying  are 
relatively  constant,  with  minor  fluctuations.  However,  if  the  Voltaire  impact  occurs 
when  Deimos  is  near  periapsis,  a  spike  in  Mars  impacts  and  a  corresponding  drop  in 
particles  escaping  the  system  are  noted.  This  is,  in  fact,  due  to  the  difference  in  true 
longitude  between  Phobos  and  Deimos  at  the  time  of  ejecta  launch;  depending  on  the 
relative  conjunctive  geometry,  the  mass  flux  from  Deimos  to  Phobos  can  be  up  to 
500%  higher. 

It  is  surprising  to  note  that  there  is  more  mass  flux  from  Deimos  to  Phobos  as 
opposed  to  Mars;  intuitively,  one  would  expect  that  most  mass  ejected  from  Deimos 
would  either  reaccrete  or  spiral  down  to  Mars.  Because  greater  mass  is  released  at 
lower  ejection  velocities  [O’Keefe  and  Ahrens,  1985],  these  results  suggest  impacts 
on  Deimos  may  have  an  effect  on  Phobos'  geology;  we  shall  attempt  to  estimate  the 
magnitude  of  that  effect  in  Section  2.5. 

For  impacts  with  Mars,  Phobos  or  Deimos,  impact  velocity  is  calculated  with 
reference  to  the  Planet-Centered  Planet  Fixed  frame  (Section  2.3. 1.3)  in  question.  The 
variation  of  impact  speeds  at  Phobos  and  Deimos  is  charted  across  28  GCs  (Figure 
2.4).  For  Phobos,  regardless  of  orbital  location  at  the  time  of  ejecta  release,  impact 
velocity  scales  linearly  with  particle  ejection  velocity.  Faster  particles  impact  with 
higher  speeds,  in  some  cases  up  to  4  km/s  (though  still  not  as  high  as  ~20  km/s 


18 


expected  for  heliocentric  impactors).  On  Deimos,  however,  impacts  above  1  km/s  are 
relatively  rare.  Almost  no  high-velocity  (>1  km/s)  impacts  are  noted  from  near- 
apoapsis  positions.  81%  of  impacts  are  clustered  in  the  0.4-0.8  km/s  region,  implying 
that  low-velocity  reaccretions  to  Deimos  are  relatively  common. 

Next,  the  relationship  between  impact  velocities  and  the  time  to  impact  is 
examined.  No  significant  acceleration  effect  with  time  or  particular  links  to  orbital 
geometry  are  noted  (Figure  2.5).  However,  the  contrast  between  the  two  bodies  is 
again  evident.  On  Phobos,  the  majority  of  impacts  occur  in  less  than  100  years; 
subsequent  impacts  become  less  frequent  as  time  increases.  This  implies  that  the  500- 
year  timescale  selected  is  adequate  to  capture  the  majority  of  Deimos-to-Phobos 
material  transfer.  On  Deimos,  however,  impacts  continue  to  build,  with  the  flux  of 
impacts  remaining  relatively  constant  even  at  the  end  of  the  500-year  timescale. 
Therefore,  it  seems  that  while  the  majority  of  mass  transfer  to  Phobos  occurs  early  on, 
reaccretions  to  Deimos  likely  continue  out  to  the  lO"^  year  timescale  hypothesized  by 
Soter  (1971).  This  also  makes  it  likely  that  a  large  number  of  the  particles  still  flying 
at  the  end  of  the  500-year  simulation  will  end  up  reaccreting  to  Deimos.  As  a 
consequence,  it  is  unlikely  that  increasing  the  computation  time  will  qualitatively 
change  our  conclusions. 

The  flight  path  angle  (FPA)  is  the  angle  between  the  incoming  velocity  vector 
and  the  position  vector  defined  by  the  surface  of  the  planet,  and  can  be  calculated  as 
[Curtis,  2013]: 


tan  Y  =  - 


2 


e  sin  O 
1+e cos  O 


(2.15) 
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where  e  is  the  eccentricity  of  the  impact  trajectory  and  (I>  is  the  true  anomaly. 

The  relationship  between  impact  velocity  and  FPA  reinforces  the  rarity  of 
high-speed  reaccretion  events  on  Deimos  (Figure  2.6).  On  Phobos,  the  frequency 
distribution  of  impactor  velocities  and  flight  path  angles  suggests  that  impacts  created 
by  Deimos  ejecta  can  vary  from  oblique,  classically  secondary  impacts  to  direct 
cratering  events. 

Finally,  investigating  the  likelihood  of  continuing  collisions  beyond  the 
chosen  500-year  timeframe  tests  the  fidelity  of  these  results.  Figure  2.7  shows  time 
curves  for  particles  still  flying  and  particles  impacting  Phobos  or  Deimos.  As 
expected,  Phobos  impacts  taper  off  with  time,  and  the  total  number  of  Phobos  impacts 
(y)  fits  well  (R^  >  0.995)  to  a  logarithmic  distribution  defined  by  y  =  178.4  ln(t/to), 
where  t  is  time  in  years  and  the  time  constant  to  is  -10.46  years.  Deimos  impacts 
continue  to  increase  and  fit  well  (R^  >  0.995)  to  a  distribution  defined  by  y  = 
4691n(t/to),  where  the  time  constant  to  is  -47.9  years.  Interestingly,  the  ratio  between 
the  time  constants  for  Phobos  and  Deimos  are  similar  to  the  ratio  of  their  orbital 
periods. 

Reaccretions  to  Deimos  are  therefore  expected  to  continue,  but  for  how  long? 

The  shape  of  the  graph  for  particles  still  in  flight  is  in  a  logarithmic  decrease;  when 

extrapolated  (Figure  2.7,  far  right)  it  takes  approximately  10,000  years  for  the  number 

of  particles  still  in  flight  to  drop  below  10%  of  the  total  number  of  particles 

generated.  This  result  agrees  well  with  predictions  made  by  Soter  (1971).  However,  a 

word  of  caution  is  appropriate  here.  Soter  and  this  study  both  neglect  effects  of  solar 
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radiation  pressure.  While  a  second-order  effect  on  the  500-year  timescale,  it  can  play 
a  significant  role  across  longer  time  periods,  and  could  decrease  the  time  to  impact 
[Klacka,  2002],  Based  on  this,  conclusions  drawn  here  are  not  expected  to  change 
qualitatively  with  an  increase  in  propagation  time  and  corresponding  decrease  in 
number  of  particles  still  flying. 

By  analyzing  the  Martian  system  within  the  framework  of  an  analytical 
restricted  three-body  problem,  previous  work  finds  that  essentially  all  ejecta  from 
either  Phobos  or  Deimos  will  be  reaccreted  to  the  moon  of  origin  \Soter,  1971; 
Dobrovolskis  and  Bums,  1980].  Differing  results  in  this  work  suggest  that  the  four- 
body  ephemerides  formulation  is  critical  to  full  understanding  of  the  orbital 
dynamics. 

2.5  Mass  Transfer  between  Martian  Satellites 

To  investigate  the  geologic  impact  of  mass  transfer  between  Phobos  and 
Deimos,  it  is  desired  to  convolve  the  probability  distributions  from  Figure  2.3  with  an 
appropriate  mass-velocity  distribution.  Advanced  scaling  laws  developed  from 
numerical  methods  exist  [e.g.  Leinhardt  and  Stewart,  2012],  but  given  the 
uncertainties  associated  with  several  key  parameters,  a  more  transparent  and  simpler 
approach  is  preferred.  For  gravity-dominated  cratering  the  volume  ejected  faster  than 
a  given  velocity  is  \Holsapple,  1993;  Housen  and  Holsapple,  2011]: 
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where  Vgj  is  the  ejection  velocity  (from  Figure  2.3)  and  R  is  the  final  radius  of 
Voltaire. 

A  sand-like  surface  is  well  represented  by  v  =  1.2  [Melosh,  1989]. 
Experimental  results  that  determine  mass-velocity  distributions  for  impacts  into 
granular  targets  find  Cgj  =  0.25  [Hermalyn  and  Schultz,  2013].  Their  results  correlate 
well  to  the  literature;  Andrews  (1975),  Cintala  et  al  (1999)  and  Stoffler  et  al.  (1975) 
find  Cgj  values  between  0.25  and  0.36.  We  adopt  Cgj  =  0.3  and  v  =  1.2.  For  Phobos' 
density  p  =1.9  g/cc  [Avanesov  et  al.,  1989;  Rosenblatt  et  al.,  2008;  Schmedemann  et 
al.,  2014]  is  used.  From  (2.16),  the  mass  ejected  within  each  velocity  bin  is  calculated 
(Figure  2.8).  In  total,  3x10^  kg  is  ejected  from  Deimos  between  400-2000  m/s. 

Using  Equation  2.16  and  the  acceleration  due  to  gravity  of  Deimos  (0.003 
ms^),  the  total  mass  excavated  faster  than  escape  velocity  is  6.1x10^'  kg;  0.5%  of  this 
total  is  therefore  ejected  in  the  400-2000  m/s  velocity  range.  Of  this  0.5%,  what 
percentage  reaches  Phobos?  From  integrating  the  mass  delivered  per  velocity  bin 
(Figure  2.8)  and  averaging  it  across  impacts  at  all  GCs  (Figure  2.3),  approximately 
3.5x10^  kg  impacts  Phobos.  This  is  12%  of  the  total  mass  released  in  the  400-2000 
m/s  range.  This  is  also  21%  of  the  mass  not  still  in  orbit  at  the  end  of  the  simulation 
(“still  flying”.  Figure  2.3).  If  all  this  mass  were  to  ultimately  impact  Phobos,  the  total 

o 

impacting  mas  would  be  6x10  kg,  the  likely  maximum  value.  Therefore,  12-21%  of 
the  mass  that  can  reach  Phobos  does  end  up  on  Phobos  on  a  10"^-year  timescale. 

The  same  analysis  yields  7.2x10^  kg  impacting  Deimos,  which  is  25%  of  the 
total  mass  released  in  the  400-2000  m/s  range  and  42%  of  the  mass  in  this  range  not 
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still  flying.  This  yields  an  upper  bound  of  1.2x10^  kg  impacting  Deimos.  Ejecta 
launched  at  lower  velocities  than  400  m/s  cannot  reach  Phobos  and  will  mostly  re¬ 
impact  Deimos.  Mars  only  receives  5-9%  of  ejecta  in  the  400-2000  m/s  range.  The 
remaining  ejecta  escapes  the  Mars  system  into  heliocentric  space. 

2.5.1.  Comparison  with  Background  Flux 

We  seek  to  place  the  sesquinary  mass  flux  into  perspective  by  comparing  it  to 
the  estimated  background  mass  flux  from  meteoroidal  impacts.  The  exact  flux  of 
small  meteoritic  bodies  at  Mars  orbit  is  not  known,  so  a  formulation  dependent  on  a 
characteristic  timescale  is  derived  below,  which  eliminates  the  Mars  mass  flux 
quantity.  Brown  et  al  (2002)  use  data  from  geostationary  satellites  around  the  Earth  to 
estimate  a  power  law  relationship  between  the  number  of  objects  colliding  with  the 
Earth  per  year  (AO  with  diameters  of  at  least  D,  of  the  form: 

log  IV  =  Co  —  dologD  (2.17) 

where  Cq  =  1.568  ±  0.03,  do  =  2.70  ±  0.08.  Assuming  the  same  power  law 
distribution  for  Martian  system  bodies,  for  N  =  Npionet  this  can  be  reformulated  as: 

^planeti.'^  ^  ~  ^planet ^  (  2.18  ) 

where  Cpianet  =  10^^“.  As  will  be  shown,  the  value  of  this  constant  does  not 
matter  for  the  analysis  in  this  work.  Assuming  spherical  impactors  with  diameter  D 
and  density  Pj,  the  incremental  number  of  impactors  dN  per  year  results  in  a  mass 
flux  increment  of: 

dM  =^PinD^dN  (2.19) 
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kg  per  year.  Differentiating  Equation  2.18,  substituting  into  Equation  2.19  and 
integrating  to  D^^ax  =  D,  the  mass  accumulated  by  a  generic  planetary  body  in  kg/yr 
from  asteroidal  flux  is: 


^0  ^Pi 
—  .  6  . 


£)3-do 


( 2.20 ) 


^planetiA  ^  ~  ^planet 

where  do  <3.  Zahnle  et  al  (2003)  derive  a  relationship  for  the  impact  rate  of  a 
satellite  compared  to  its  planet.  Applied  to  the  Mars  system,  this  is: 

^Pho(A  ^  ^marsfpho  (  2.21  ) 

where  D  is  the  diameter  of  the  largest  impactor  incident  to  Phobos  and: 

fpho  =  (  2.22  ) 

^mars^Pho 

-7  -f\ 

where  a  is  the  distance  from  Mars.  The  resulting  values  are  5x10'  and  4x10' 
for  Deimos  and  Phobos  respectively.  Combining  Equations  2.18  and  2.21  for  Phobos: 

( 2.23  ) 


n  _  f Pho  jijst 

^Pho  ~  „-do  ‘^mars 
^st 


where  IV^ars  signifies  >  D^^).  The  assumption  here  is  that  Stickney 

is  the  largest  impact  to  have  occurred  on  Phobos  over  its  history  and  the  diameter  of 
the  Stickney-forming  impactor  is  Substituting  Equation  2.23  into  Equation  2.18, 
and  applying  2.18  to  N^ars'- 


Mno(d  >  Ds,)  =  C^arsfpi,o  [^]  [^]  d’,'"”  (  2  24  ) 

Next,  we  define  a  characteristic  timescale  t  defined  such  that  t  years  elapse 
between  Voltaire-size  collisions  on  Deimos.  Erom  (2.21)  and  the  definition  of  N: 

T  =  =  l/ifoeiK^drs)  (  2.25  ) 
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where  ^mars(.d  >  £>voitaire)  foi"  Deimos.  Expanding  this  according 

to  Equation  2.18: 

T  =  )  (  2.26  ) 

Across  T  years,  combining  Equations  2.24  and  2.26,  the  poorly  known  (and, 
on  long  timescales,  time- variable)  annual  mass  flux  delivered  to  Mars  represented  by 
Cmars  cancels.  The  total  mass  accreted  by  Phobos  due  to  background  impacts  on  a 
Voltaire  timescale  is: 


'^PhOA 


T  MphQ  — 


f  P  ho 
f  Dei 


r  do  1  [zEEil 
b-dol  [  6  J  D-^0 


( 2.27 ) 


Equation  2.27  is  only  dependent  on  the  diameter  of  the  impactors  and  the 
slope  of  the  size-frequency  distribution.  Eor  Phobos,  the  largest  crater  is  Stickney, 
with  a  170-m  likely  impactor  size  [Asphaug  and  Melosh,  1993].  For  Deimos,  the 
Voltaire  impactor  diameter  is  estimated  from  gravity-dominated  scaling  relations, 
rearranged  from  Cintala  and  Grieve  (1998)  and  Schmidt  and  Housen  (1987)  and 
similar  to  Zahnle  et  al  (2003): 

.  1  ,  1.2821 


=  (  0.862Z), 


( 2.28  ) 


where  the  subscript  i  denotes  the  impactor  that  created  Voltaire,  subscript  d 
denotes  Deimos,  units  are  CGS  and  D,  is  the  diameter  of  the  Voltaire  transient  crater, 
taken  to  be  1.95  km.  Asphaug  and  Melosh  (1993)  assume  an  impact  of  3  km/s  for  the 
Stickney  impact;  assuming  the  same  impact  velocity  and  impactor  density,  (2.28) 
yields  a  Voltaire  impactor  diameter  of  25  m.  From  Equation  2.27,  for  an  asteroid-type 


impact  (pj=2.6  g/cc;  Barnhart  and  Nimmo  (201 1)),  the  mass  accreted  by  Phobos  from 
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solar  system  impactors  between  Voltaire-scale  impacts  is  approximately  2.9x10^  kg; 
it  would  be  less  if  the  impactor  were  assumed  to  have  the  lower  characteristic  density 
of  Phobos.  Comparing  this  to  the  3.5x10^  kg  that  impacts  Phobos  during  every 
Voltaire  impact,  the  fraction  of  Deimos  material  delivered  to  Phobos  (F)  represents, 
on  average,  0.12  of  the  material  accreted  to  Phobos  from  direct  solar  system  impactor 
flux. 

Using  a  similar  derivation  for  Deimos,  it  is  found  that  only  2x10^  kg  is 
accreted  to  Deimos  by  solar  system  impactors  between  Voltaire-size  events.  Deimos 
receives  less  material  than  Phobos  because  the  focusing  factor  (Equation  2.21)  is 
smaller.  This  flux  is  exceeded  greatly  by  the  reaccreted  mass  ejected  during  a 
Voltaire-size  impact;  for  ejecta  with  velocities  400-2000  m/s,  F  =  3.6,  where  F  is  the 
fraction  of  Deimos  material  delivered  to  Phobos.  However  this  velocity  range  is  only 
0.5%  of  the  total  mass  thrown  into  orbit.  Ejecta  with  velocity  <  400  m/s  has 
insufficient  velocity  to  reach  Phobos  (or  Mars)  and  will  likely  reaccrete  to  Deimos. 
Thus,  the  true  value  of  F  across  all  ejecta  is  likely  -700  for  Deimos.  In  other  words, 
on  Deimos,  mass  reaccreted  during  a  Voltaire-size  impact  greatly  exceeds  mass 
naturally  accreted  from  solar  system  impactors.  However,  sesquinaries  from  a  large 
impact  on  Deimos  provide  only  a  minor  contribution  to  the  flux  at  Phobos. 

As  a  reality  check,  the  applicability  of  Equation  2.17  to  Mars  system  impacts 
is  calculated.  The  frequency  of  a  Voltaire-sized  impact  is  t  =  MNdei  years,  or  134  Ma. 
The  issue  of  the  apparent  age  of  Voltaire  is  discussed  further  below.  Using  Equation 
2.26,  the  ratio  of  the  timescale  for  a  Voltaire-forming  impact  on  Deimos  to  a 
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Stickney-forming  impact  on  Phobos  is  0.049.  Therefore,  the  frequency  of  a  Stickney- 
sized  impact  is  approximately  every  2.7  Ga,  which  seems  quite  reasonable. 

As  an  additional  check,  the  mass  flux  to  Mars  is  estimated  from  studies  in  the 
literature.  Chappaz  et  al  (2011)  find  that  the  mass  flux  from  Mars  to  Phobos  can  be 
estimated  at  0.25  jUg/m^/yr  or  0.2  kg/yr  across  the  inner  moon’s  surface  area  assuming 
a  mixing  depth  of  0.5  m.  The  mass  flux  of  solar  system  projectiles  to  Phobos  is 
numerically  found  to  be  k  =  40-2400  times  greater  than  the  flux  from  Mars  ejecta, 
with  a  Monte  Carlo  preferred  value  of  k  =  195  [Ramsley  and  Head,  2013].  The 
preferred  value  yields  a  mass  flux  at  Phobos  of  38.5  kg/yr  from  asteroids,  comets  and 
meteoroids.  At  this  rate,  across  134  Ma,  Phobos  accumulates  5.2x10^  kg  from  direct 
impacts,  within  a  factor  of  two  of  the  2.9x10^  kg  derived  above. 

While  independent  of  Mars’  meteoroidal  flux,  the  results  are  admittedly 
susceptible  to  the  assumed  size  of  the  Voltaire  impactor.  The  size  of  Stickney  implies 
that  the  gravity  regime  approximation  is  likely  appropriate.  However  the  short 
timescale  derived  for  the  Voltaire  impact  implies  that  it  may  be  deeper  in  the  strength 
regime  than  initially  assumed.  If  true,  the  Voltaire  impactor  would  have  to  be  larger 
than  25  m,  which  would  reduce  quantitative  estimates  of  F,  the  fraction  of  Deimos' 
mass  flux  with  respect  to  the  solar  system  impactor  mass  flux. 

Assuming  that  the  impactor  was  incident  at  a  45°  angle  [Melosh,  1989; 
Holsapple,  1993],  the  diameter  of  an  impactor  for  strength-dominated  cratering  may 
be  estimated  by  rearranging  Equation  1.17  from  Zahnle  et  al  (2008)  to  yield: 
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D„„,  =  1.027  (2.29) 

where  Y  represents  the  dynamic  strength  of  the  body;  other  variables  are  as  in 
Equation  2.28.  Melosh  (1989j  uses  a  value  of  F  =  2  MPa,  the  observed  yield  stress  at 
crater  collapse,  to  estimate  the  gravity/strength  transition  on  the  Earth  and  the  Moon. 
Adopting  this  value  in  Equation  2.29,  the  diameter  of  the  Voltaire  impactor  is  88  m. 
However,  the  value  of  Y  is  uncertain;  for  the  low-density  Deimos,  it  is  unlikely  that 
the  yield  strength  is  as  high  as  2  MPa.  Eor  an  order  of  magnitude  change  in  F  ranging 
from  0.2  to  20  MPa,  the  diameter  of  the  impactor  varies  from  40  to  190  m. 

Is  strength  or  gravity  scaling  more  appropriate  for  modeling  the  Voltaire 
impact?  Gravity  can  be  a  factor  on  solar  system  bodies  as  small  as  400  m  [Love  and 
Ahrens,  1996].  Phobos  and  Deimos  have  similar  compositions,  bulk  densities  and 
accelerations  due  to  gravity  [Davis  et  al,  1981;  Szeto,  1983],  and  several  arguments 
for  modeling  cratering  on  Phobos  in  the  gravity-dominated  regime  are  detailed  in 
Asphaug  and  Melosh  (1993)  and  Asphaug  et  al  (2015a).  Finally,  while  some  authors 
have  used  the  wide  distribution  of  ejecta  on  Deimos  to  surmise  that  strength-scaling 
may  be  appropriate  for  Deimos  [Lee  et  al,  1986],  we  suggest  that  this  global 
distribution  could  instead  be  a  function  of  the  large  amount  of  mass  reaccreted  over 
500yr  timescales  (Section  2.3).  However,  even  if  the  Voltaire  impact  is  in  the  strength 
regime,  this  would  not  change  our  qualitative  results,  i.e.,  1 )  that  the  sesquinary  mass 
transfer  is  a  relatively  small  fraction  of  meteoroidal  impacts  to  Phobos,  and  2)  that  the 
mass  reaccreted  to  Deimos  from  large  Deimos  impacts  exceeds  the  meteoroidal  mass 
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flux  to  Deimos.  For  instance,  increasing  the  diameter  of  the  Voltaire  impactor  to  88 
m,  the  qualitative  results  become  F  =  0.004  for  Phobos  and  F  =24  for  Deimos  across 
all  ejecta  velocities.  Therefore,  uncertainty  on  where  the  Voltaire  impact  lies  in  the 
gravity/strength  regime  does  not  affect  our  qualitative  conclusions. 

2.6.  Discussion 

2.6.1.  Importance  of  Inter-Moon  Mass  Transfer  and 
Reaccretion 

The  central  result  of  this  work  concerns  the  relative  mass  transfer  during 
impacts  on  Deimos.  It  has  been  found  that  a  Voltaire-sized  impact  on  Deimos  does 
transfer  mass  from  Deimos  to  Phobos,  with  sufficient  velocity  to  create  primary  crater 
morphology,  discussed  further  below.  When  viewed  with  reference  to  the  solar 
system  impactor  flux,  the  sesquinary  mass  transfer  is  not  significant,  and  is  likely  in 
the  10%  range.  However,  compared  to  the  25  ppm  of  Mars  material  in  Phobos 
regolith  estimated  by  Chappaz  et  al  (2011),  one  out  of  ten  particles  originating  from 
Deimos  is  a  large  number  and  is  of  interest  to  Phobos  lander  mission  concepts  [e.g. 
Udrea  et  al.,  2015,  2016]. 

While  the  placement  of  Voltaire  within  the  strength-gravity  domain  can  cause 
some  uncertainty  in  impactor  size,  it  has  been  shown  that  the  total  ejecta  mass 
reaccreted  to  Deimos  is  likely  to  greatly  exceed  the  background  mass  flux  to  Deimos 
(F  >  20).  A  Voltaire- sized  impact  could  therefore  feasibly  resurface  large  parts  of  the 
moon,  erasing  the  previous  geological  record.  Dating  the  surface  of  Deimos  may  be 
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more  challenging  than  previously  suspected;  the  surface  age  may  better  represent  the 
age  of  either  Voltaire  or  the  similarly  sized  Swift  crater.  Further,  an  11 -km  concavity 
on  the  southern  end  of  Deimos  is  hypothesized  to  be  a  possible  impact  scar  from  an 
ancient,  very  large  collision  [Lee  et  al,  1986;  Thomas,  1993;  Thomas  et  al,  1996].  If 
true,  this  would  have  resulted  in  the  transfer  of  significant  sesquinary  mass  transfer  to 
Phobos,  and  a  complete  resurfacing  of  Deimos’ s  surface. 

2.6.2  The  Spectral  Dichotomy  of  Phobos 

Phobos  exhibits  two  distinct  spectral  units,  one  of  “redder”  origin  and  one  of 
“bluer”  origin,  a  distinction  that  likely  stems  from  a  compositional  difference  [Rivkin 
et  al,  2002;  Lee,  2009].  It  has  been  proposed  that  the  red  unit  is  a  wide-spread 
shallow  layer  superimposed  on  a  blue  base  that  is  perhaps  more  representative  of 
Phobos  composition  at  depth  [Murchie  and  Erard,  1996;  Murchie  et  ah,  2008]. 
Possible  mechanisms  for  the  superimposition  of  red  material  were  briefly  outlined  in 
Section  2.1  [Britt  and  Pieters,  1988].  These  include  the  hypothesis  that  the  spectrally 
red  “veneer”  may  be  ejecta  accreted  from  Deimos  rather  than  Mars,  as  suggested  by, 
e.g..  Smith  et  al  (2015). 

Presented  results  for  the  distribution  of  low-velocity,  oblique-angle  impacts 

(Figure  2.6)  support  the  existence  of  trajectories  that  could,  in  theory,  deposit  a 

“veneer”  of  red  Deimos  material  across  Phobos’ s  surface.  However,  the  inter-moon 

mass  flux  is  relatively  small  compared  to  the  solar  system  impactor  flux,  which  likely 

has  a  greater  effect  on  the  global  surface  geology,  particularly  in  the  lOO-t  Ma  since 

the  last  Voltaire- sized  impact.  Therefore  it  is  believed  to  be  unlikely  that  the  red 
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veneer  of  Phobos  is  of  Deimos  origin.  Recent  spectral  analysis  by  Thomas  et  al 
(2011)  supports  this  finding  with  evidence  of  subsequent  impacts  penetrating  the  blue 
unit  near  Stickney  to  reveal  redder  material.  This  observation  casts  doubt  on  the  idea 
that  the  blue  unit  may  be  representative  of  Phobos  at  depth  [Basilevsky  et  al,  2014]. 
Ultimately,  sample  return  from  both  Phobos  and  Deimos  will  conclusively  establish 
which  spectral  unit  is  representative  of  depth;  results  from  this  work  suggest  that  the 
surface  unit  is  unlikely  to  originate  from  Deimos. 

2.6.3  Primary  Versus  Secondary  Impact  Morphology 

Due  to  the  dearth  of  classic  secondary  impact  features  such  as  radial  crater 
chains  or  herringbones,  it  has  been  concluded  that  few,  if  any  Phobos  craters  are 
secondary  in  origin  [Thomas,  1979;  Thomas  and  Veverka,  1980a].  A  similar 
conclusion  was  reached  for  Deimos.  A  limit  of  10  m/s  on  maximum  re-impact 
velocity  was  proposed  by  Thomas  (1998)  and  continues  to  be  used  in  analysis  of 
Phobos'  geology  [Murray  et  al,  2006;  Schmedemann  et  al,  2014;  Smith  et  al,  2015]. 
However,  every  sesquinary  impact  studied  here  occurs  at  speeds  above  100  m/s. 
Though  still  orders  of  magnitude  below  heliocentric  impactor  velocities.  Figure  2.6 
shows  several  high-velocity,  high-FPA  particles  that  could  create  craters 
indistinguishable  from  primary  impact  craters  on  Deimos,  and  particularly  on  Phobos. 

On  Deimos,  over  80%  of  impacts  cluster  in  the  0.6  ±  0.2  km/s  region, 
implying  that  the  majority  of  reaccretions  are  low  velocity  (subsonic).  However  for 
Phobos  sesquinary  particles  can  arrive  with  either  subsonic  or  supersonic  impact 
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velocities.  One  would  therefore  predict  a  wide  range  in  the  resulting  crater 
morphology;  a  comprehensive  image  survey  of  Phobos  may  be  able  to  distinguish 
between  the  different  crater  morphologies.  Several  low-velocity,  low-FPA  particles 
incident  to  Phobos  have  been  found  here,  which  could  create  classic  oblique  or 
chained  secondary  impact  morphology;  a  Phobos  image  survey  by  Smith  etal  (2015) 
finds  several  craters  and  deposits  likely  originating  from  such  low-velocity  impacts. 

The  escape  speed  from  Mars  at  Phobos’  orbit  is  3  km/s,  which  is  exceeded  by 
the  fastest  sesquinary  impacts  (Figure  2.4).  These  impacts  can  create  ejecta  of  their 
own  that  may  subsequently  be  lost  from  the  Mars  system.  Mass  loss  from  Phobos  is 
one  possible  explanation  for  why  outlines  of  ejecta  blankets  are  not  conspicuous  on 
Phobos  [Lee  et  al,  1986].  Additional  simulations  of  Phobos-centered  ejecta  dynamics 
are  presented  in  Chapter  3  and  confirm  this  hypothesis. 

2.6.4  Surface  smoothness,  Brightness  and  Impact  gardening 

Very  little  ejecta  escapes  large  bodies  (e.g.  Earth),  with  the  majority 
redeposited  locally  as  a  continuous  ejecta  blanket.  Most  ejecta  escape  from  very  small 
bodies  (e.g.  Phoebe  [Burns  et  al,  1996]),  never  to  be  reaccreted.  Ejecta  dynamics  on 
Deimos  present  an  interesting  bridge  between  these  two  regimes,  with  sesquinary 
effects  appearing  to  be  important.  Ejecta  escapes  but  is  then  reaccreted  on  a  timescale 
of  up  to  hundreds  of  years,  resulting  in  a  global,  near-isotropic  redistribution  of 
sesquinary  ejecta. 
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For  a  Voltaire- sized  impact,  almost  all  of  the  ejecta  material  launched  at 
velocities  <  400  m/s  will  ultimately  reaccrete  to  Deimos  [Soter,  1972;  Thomas  etal, 
1996].  From  equation  (2.16)  above  this  represents  6.1xl0"  kg,  or  about  0.5m 
thickness  of  material  distributed  evenly  over  the  entire  surface.  We  estimate  that  the 
Voltaire  and  Swift  impacts  together  could  have  added  on  the  order  of  1  m  of  fresh 
regolith  or  crater  fill  to  Deimos.  It  has  been  suggested  that  the  smoother  and  brighter 
surface  of  Deimos  is  due  to  crater  fill  of  5-7  m  depth  [Veverka,  1978;  Thomas,  1979, 
1993;  Thomas  et  al,  1996].  Given  our  estimate,  it  is  suggested  that  reaccreting 
sesquinary  mass  provides  at  least  a  partial  explanation  for  the  origin  of  this  crater  fill 
material. 

The  impact  velocity  distributions  (Figure  2.4)  also  suggest  a  tertiary  ejecta 
effect.  Though  most  sesquinary  impacts  to  Deimos  are  relatively  low  speed  (<1 
km/s),  these  are  still  significantly  higher  than  the  escape  velocity  of  Deimos  (5.5  m/s. 
Table  2.1)  and  could  potentially  launch  additional  ejecta  in  their  turn  (see  below). 
Similarly,  ejecta  from  Deimos  can  impact  Phobos  with  enough  mass  and  speed  to 
create  craters  and  excavate  Phobos  mass,  which  could  then  enter  Mars  orbit. 

On  Deimos,  sesquinary  impacts  represent  a  large  mass  flux  relative  to  the 
background  flux  and  have  ~km/s  impact  velocities  (Figure  2.2).  The  result  is  likely  to 
be  the  production  of  further  suborbital,  ballistically  emplaced  ejecta.  This  mechanism 
could  contribute  to  the  smooth  appearance  of  the  Deimos  surface.  Energy  and 
momentum  transfer  from  reaccretion  impacts  might  even  set  off  downslope 
movement  noted  on  Deimos  [Thomas  and  Veverka,  1980b],  though  admittedly  the 
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efficiency  of  the  impacts  at  initiating  this  process  is  unknown  and  should  be 
constrained  in  the  future. 

The  possibility  of  sesquinary  impact  gardening  makes  the  evolution  of 
regolith  in  between  major  collision  events  a  complex  process  on  both  Martian  moons. 
Determining  the  exact  nature  of  Deimos  material  mixed  with  Phobos  regolith  is  one 
of  the  primary  science  objectives  of  lander  concepts  in  development  for  Phobos 
\Udrea  et  a/.,  2015,  2016].  The  methods  applied  in  this  work  are  further  applicable  to 
understanding  regolith  development  and  dust  belts  on  small  bodies  within  planetary 
gravitational  wells,  such  as  the  Saturnian  moons  of  Tethys,  Calypso  and  Telesto. 

2.7  Figures 


Figure  2.1.  Illustration  of  the  relationship  between  the  Planet-Centered  Inertial  (PCI)  and 
Planet-Centered  Planet-Fixed  (PCPF)  frames. 
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Figure  2.2.  Illustration  of  the  relationship  between  the  Planet-Centered  Planet  Fixed  (PCPF) 
frame  and  the  South-East-Zenith  (SEZ)  frames.  The  fundamental  plane  is  highlighted  (blue). 


Fate  of  ejecta:  Function  of  orbital  position 
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Figure  2.3.  The  fate  of  Deimos  ejecta  from  Voltaire  as  a  function  of  (Left)  Orbital  Geometry  and 
(Right)  Ejection  Velocity  from  Voltaire.  GC-1  is  near-apoapsis;  GC-14  is  near-periapsis.  Note 
that  the  difference  in  the  Deimos-Phobos  angle  is  the  important  quantity  with  regard  to  ejecta 
fate  (see  text). 
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Figure  2.4.  (Left)  Impact  Velocity  at  Phobos  and  (Right)  Impact  Velocity  at  Deimos  versus 
particle  ejection  velocity  from  Voltaire.  While  Phobos  exhibits  several  high-velocity  impacts, 
high-velocity  impacts  at  Deimos  are  relatively  rare  and  are  primarily  clustered  below  1  km/s 
(compare  to  -20  km/s  heliocentric  impactor  velocity).  The  discontinuity  at  700  m/s  in  both 
graphs  is  due  to  an  increase  in  impacts  to  Mars  at  that  velocity  range  for  certain  Deimos-Phobos 
orbital  geometries  (Figure  2.3). 
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Figure  2.5.  (Left)  Impact  velocity  at  Phobos  and  (Right)  Impact  velocity  at  Deimos  versus 
duration  of  particle  flight.  While  impacts  to  Phobos  are  frequent  in  the  first  100  years  post 
impact,  they  begin  to  taper  off  toward  the  end  of  the  examined  duration;  impacts  to  Deimos,  on 
the  other  hand,  continue  at  a  relatively  constant  pace,  implying  that  while  500  years  captures  the 
bulk  of  Phobos  mass  transfer,  reaccretions  to  Deimos  will  likely  continue  to  the  10"^  year 
timescale  [Soter,  1971;  Davis  etaL,  1981]. 
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Phobos  impact  dynamics:  Angle  vs  velocity 


Deimos  impact  dynamics:  Angle  vs  velocity 
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Figure  2.6.  (Left)  Impact  flight  path  angle  (FPA)  at  Phobos  and  (Right)  Impact  FPA  at  Deimos 
versus  the  impact  velocities  on  the  respective  bodies.  Apart  from  reinforcing  the  fact  that  Deimos 
impacts  are  primarily  low-velocity,  we  also  see  a  wide  distribution  in  flight  path  angles.  There 
are  several  low-velocity,  low-FPA  impacts  that  should  create  oblique  or  secondary  crater 
morphology,  and  several  high-velocity,  high-FPA  impacts  that  will  exhibit  direct  or  primary 
crater  morphology. 


Figure  2.7.  Number  of  (Left)  particles  still  flying  and  (middle)  impacts  to  Phobos  and  Deimos 
with  time.  The  curves  are  well  behaved,  with  no  unexpected  jumps.  Impacts  to  Phobos  can  be 
seen  to  be  tapering  off,  while  reaccretions  to  Deimos  continue  to  rise.  These  are  expected  to 
continue  until  no  more  particles  are  still  flying.  (Right)  Extrapolation  of  the  logarithmic  decrease 
in  the  left  plot.  The  decrease  fits  well  (R^  >  0.995)  to  a  logarithmic  distribution  defined  by  y  =  - 
5631n(t/to),  where  to  --25,000  years.  Using  this  curve,  it  takes  10,000  years  for  the  particles  still 
flying  to  drop  below  500,  agreeing  with  Soter  (1971). 
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Ejection  velocity  (m/s) 


Top  of  velocity  bin  (100  m/s  intervals) 


Figure  2.8.  (Left)  Total  mass  faster  than  ejection  velocities  (Equation  2.16).  (Right)  Mass  ejected 
per  velocity  bin  in  tens  of  millions  of  kg.  Velocity  bins  correspond  to  Figure  2.3. 
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Chapter  3 

Sesquinaries:  A  Case  Study  of  Phobos 


This  chapter  is  a  modified  reprint  ofNayak  and  Asphaug  (2016),  Sesquinary  Catenae 
on  the  Martian  Satellite  Phobos  from  Reaccretion  of  Escaping  Ejecta,  Nature 
Communications,  DOI:  10.1038/ncommsl2591. 

3.1  Abstract 

The  Martian  satellite  Phobos  is  crisscrossed  by  linear  grooves  and  crater 
chains  whose  origin  is  unexplained.  Anomalous  grooves  are  relatively  young  and 
crosscut  tidally  predicted  stress  fields  as  Phobos  spirals  toward  Mars.  In  this  chapter, 
we  report  strong  correspondence  between  these  anomalous  features  and  reaccretion 
patterns  of  sesquinary  ejecta  from  impacts  on  Phobos.  Escaping  ejecta  persistently 
imprint  Phobos  with  linear,  low-velocity  crater  chains  (catenae)  that  match  the 
geometry  and  morphology  of  prominent  features  that  do  not  fit  the  tidal  model.  This 
work  proves  these  cannot  be  older  than  Phobos’  current  orbit  inside  Mars’  Roche 
limit.  Distinctive  reimpact  patterns  allow  sesquinary  craters  to  be  traced  back  to  their 
source,  for  the  first  time  across  any  planetary  body,  creating  a  novel  way  to  probe 
planetary  surface  characteristics.  For  example,  it  is  shown  that  catena-producing 
craters  likely  formed  in  the  gravity  regime,  providing  constraints  on  the  ejecta 
velocity  field  and  knowledge  of  source  crater  material  properties. 


39 


3.2  Introduction 


Ejecta  escaping  from  an  impact  on  a  natural  satellite  often  goes  into  orbit 
about  the  primary  and  can  reimpact  the  satellite  or  its  companions  after  an  extended 
period.  These  ‘sesquinary’  impacts  [Zahnle  et  al,  2008]  are  slower  than  primary 
impacts,  but  faster  than  the  satellite  escape  velocity.  Because  they  spend  time  in  orbit 
they  do  not  simply  radiate  from  their  primary  crater  like  secondary  craters; 
nonetheless  their  close  dynamical  association  with  the  satellite  can  give  them  a 
unique  geometrical  distribution.  Like  secondary  craters,  sesquinaries  are  probes  of  the 
primary  ejection  process,  but  are  also  bound  to  the  dynamics  of  the  planet-satellite 
system.  Unlike  secondaries,  to  date  no  sesquinary  crater  has  been  traced  back  to  its 
primary. 

Phobos,  the  26x22x18  km  battered  moon  of  Mars,  is  covered  in  parallel  linear 
features  whose  orientation  is,  for  the  most  part,  aligned  with  de-orbiting  tidal  stresses 
as  Phobos  spirals  closer  to  Mars.  As  the  tidal  bulge  grows,  surface  stresses  increase 
and  cause  striations  \Soter  and  Harris,  1977;  Asphaug  et  al.,  2015b].  However,  many 
of  Phobos’  linear  features  do  not  align  with  any  interpretation  of  tidal  stress,  giving 
rise  to  alternative  models  such  as  impact  fractures  related  to  the  formation  of  Stickney 
[Fujiwara,  1991],  the  largest  crater  on  Phobos,  ejecta  from  Mars  basin  formation 
[Murray  and  Heggie,  2014],  pitted  regolith  from  bouncing  boulders  [Wilson  and 
Head,  2015]  and  drainages  opening  up  into  a  fractured  suhstmte[Horstman  and 
Melosh,  1989].  Even  with  recent  improvements  incorporating  two-layer  tidal  stresses 
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from  orbital  decay  [Asphaug  et  al,  2015b],  a  subset  of  prominent  hemispherical 
crater  chains  crosscut  predicted  stress  fields  and  bear  a  closer  resemblance  to  distal 
secondary  crater  chains  on  the  Moon,  except  there  is  no  apparent  source  crater.  They 
also  resemble  the  tidal  catenae  on  Ganymede  and  Callisto  that  are  the  result  of 
cometary  disruption  inside  the  Roche  limit  of  Jupiter  [Melosh  and  Schenk,  1993].  We 
propose  these  features  on  Phobos  are  a  novel  kind  of  structure  intermediate  between 
these  two  phenomena,  which  we  call  “sesquinary  catenae”. 

We  classify  four  kinds  of  satellite  impacts.  A  primary  impact  is  by  a  bolide 
from  outside  the  planet's  sphere  of  influence;  these  are  generally  the  fastest.  A 
secondary  impact  is  by  a  fragment  ejected  from  a  primary  crater;  these  are  the  slowest 
cratering  events,  no  faster  than  a  fraction  of  the  satellite's  escape  velocity  Vesc-  Such 
craters  radiate  from  the  primary  and  form  linear  chains  and  clumps.  A  sesquinary 
impact  is  formed  by  crater  ejecta  that  escape  but  remains  in  orbit  in  the  system;  its 
impact  velocity  is  intermediate  between  Vesc  and  the  orbital  velocity  Vorb-  When  the 
satellite  is  far  from  the  planet,  sesquinaries  can  produce  primary-like  crater 
morphology.  Intermediate  between  sesquinary  and  secondary  is  the  so-called 
dosquinary,  where  the  ejecta  reimpacts  after  spending  only  a  few  orbits  aloft.  These 
can  be  thought  of  as  the  slowest  possible  sesquinaries,  but  not  quite  secondaries,  since 
the  gravitational  influence  is  primarily  that  of  the  planet.  These  are  limited  to 
satellites  that  orbit  close  to  the  planet. 

Sesquinaries/dosquinaries  from  Phobos  are  especially  interesting  for  three 
reasons:  First,  the  escape  velocity  Vesc  is  of  order  ~10  ms  ',  hundreds  of  times  slower 
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than  its  orbital  velocity  Vorb~  2.2  km  s'\  The  slowest  escaping  ejecta,  which  comprise 
the  major  mass  fraction,  cannot  stray  far  from  the  satellite.  Second,  the  escaping 
ejecta  are  subject  to  powerful  orbital  and  tidal  distortion,  since  the  current  semimajor 
axis  a~2.11  Rm^s  is  significantly  inside  the  classical  Roche  limit  Rroche~3.19  /?Mars 
[Witasse  et  al,  2014].  Third,  as  shall  be  shown,  the  timescale  for  reaccretion  is  so 
short  that  ejected  particles  can  be  reaccreted  before  they  have  time  to  disperse.  This 
leads  to  the  curious  geomorphic  phenomenon  of  sesquinary  catenae,  each  linked  to  a 
particular  source  crater  on  Phobos,  and  sets  these  features  apart  from  crater  chains  on 
bodies  beyond  deep  planetary  gravity  wells  of  planets,  such  as  asteroids  Eros 
[Buczkowski  et  al,  2008],  Gaspra  and  Steins  [Marchi  et  al,  2010].  Analytical 
formulations  have  been  used  to  study  the  dynamics  of  escaping  and  reaccreting  ejecta 
[Soter,  1971;  Bums,  1972;  Dobrovolskis  and  Burns,  1980],  but  to  precisely  predict 
sesquinary  reimpact  locations  an  ephemerides  formulation  is  required,  as  in  Chapter 
2.  Similar  to  that  approach,  multiple  orbital  configurations  around  Mars  are 
considered  to  evaluate  variations  in  Mars-Phobos  orbital  geometry  and  inferior  or 
superior  conjunctions  with  Deimos  at  the  time  of  the  primary  ejecta-producing 
impact.  Since  the  significant  percentage  of  mass  escaping  Phobos  is  ejected  slower 
than  100  ms  '  (92%,  Section  3.3),  the  component  ejected  with  velocities  from 
Vesc-llms  '  to  100  ms  '  is  tracked  in  precise  detail.  Planetocentric  latitude  and 
longitude  of  each  impact  location  are  then  extracted  from  the  evolution  model. 
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3.3  Methods 

3.3.1.  Mars  Gravity  System  Integrator 

Ejecta  are  integrated  in  an  ephemeris-centered  frame  for  an  arbitrarily  chosen 
primary  impact  epoch,  incorporating  solar  third  body  perturbations  and  a  20x20  Mars 
gravity  harmonic  model  [Murchie,  2010].  As  described  in  greater  detail  in  Chapter  2, 
our  model  accounts  for  changes  in  the  rotation  and  longitude  rate  of  Phobos, 
precession  of  the  argument  of  periapsis  and  longitude  of  the  ascending  node,  and  its 
non-spherical  shape  (triaxial  ellipsoid  with  semi-major  axes  of  13  x  11.4  x  9.1  km 
{Murchie  et  al,  2003]).  Integration  time  is  capped  at  10  years  to  permit  use  of  a  7*- 
order  Runge-Kutta  integrator  without  excessive  approximations  to  the  perturbed 
Hamiltonian  [Leimkuhler  and  Reich,  2004]. 

3.3.2.  Outbound  Ejecta  Velocity  Distribution 

For  the  ejected  material,  streamlines  ejected  by  a  primary  impact  on  Phobos 
are  modeled  using  an  axisymmetric  Z-model  formulation  [Maxwell,  1977;  Roddy, 
1977].  Though  limited  by  its  neglect  of  interactions  across  streamlines,  the  Z-model 
reasonably  approximates  several  experimentally  observed  features  [Melosh,  1989; 
Barnhart  and  Nimmo,  2011].  The  limitations  of  a  Z-model  implementation  are 
discussed  at  length  in  [Barnhart  and  Nimmo,  2011].  As  in  Chapter  2,  this  application 
is  only  concerned  with  ejecta  streamlines  that  escape  Phobos,  and  is  unaffected  by  the 
details  of  cratering  flow  beneath  the  ground  plane,  surface  material  mixing  during 
ejection  or  direct  retention  and  emplacement  of  deposits  [Nayak  et  al,  2016b]. 
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Therefore  the  Z-model  provides  a  suitable  level  of  insight  into  an  outbound  velocity 
distribution;  approximations  made  are  unlikely  to  alter  quantitative  results. 
Hydrocode  studies  of  the  Stickney  impact  (~10  km  diameter)  find  Z  =  3.5  to  suitably 
represent  cratering  flow  [Asphaug  and  Melosh,  1993],  adopted  here.  All  streamlines 
are  ejected  at  a  constant  angle  set  according  to  Equation  2.6  [Maxwell,  1977];  for  Z  = 
3.5,  this  angle  is  56.3°,  measured  from  the  horizontal  (Figure  2.1  and  2.2). 

Outbound  radial  (Vf)  and  vertical  (Vz)  ejection  velocities  vary  inversely  with 
distance  from  the  center  of  the  crater  r  as  in  Equations  2.7  and  2.8  [Maxwell,  1977], 
using  the  acceleration  due  to  gravity  g  on  Phobos  as  0.0057  ms'^.  The  cratering  flow 
is  ultimately  dependent  on  the  transient  crater  radius,  calculated  as  Ri  =  0.65  Rf,  where 
Rf  is  the  final  radius  of  the  crater  [Barnhart  and  Nimmo,  2011].  For  example,  the 
Stickney  primary  impact  is  simulated  in  this  work  using  a  final  radius  of  5.05 
km[Veverka  and  Duxbury,  1977]  and  a  transient  crater  radius  of  3.28  km.  The  number 
of  streamlines  n  is  chosen  to  yield  a  suitably  dense  distribution  of  velocities  greater 
than  the  Phobos  escape  velocity.  Setting  R^m  =  0  and  varying  i?min  <  r  <  Ri,  n 
streamlines  evenly  spaced  in  radius  are  generated  within  a  crater  of  choice. 
Streamlines  may  then  be  converted  to  axisymmetric  coordinates;  radial  and  vertical 
coordinates  are  defined  by  Equations  2.10  and  2.1 1  [Barnhart  and  Nimmo,  201 1]. 

The  Z-model  yields  a  2D  axisymmetrical  velocity  distribution;  similar  to 
Chapter  2,  this  is  then  rotated  around  the  azimuthal  direction  at  5-30°  azimuth  {P) 
intervals  to  create  a  3D  velocity  distribution.  This  three-dimensional  outbound 
velocity  distribution  is  tied  to  the  crater  being  studied,  in  the  Topocentric  Horizon 
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frame,  adapted  from  the  South-East-Zenith  (SEZ)  frame  [Vallado,  2013]  (Eigure  2.1 
and  2.2).  This  frame  is  centered  at  the  crater,  with  the  x-axis  pointing  south  and 
aligned  with  the  meridian  passing  through  the  crater.  The  x-y  plane  is  tangent  to  the 
surface  and  points  along  the  local  latitude  circle.  Completing  the  right-handed  system, 
the  z-axis  points  radially  outward  from  the  impact  site  towards  the  “local”  zenith 
[Nayak  et  al,  2016b].  Eor  use  with  the  Mars  gravity  system  integrator,  these 
coordinates  are  then  rotated  into  the  Mars  Centered  Inertial  (MCI)  frame  through  the 
Phobos  Centered  Phobos  Eixed  (PCPE)  and  Phobos  Centered  Inertial  (PCI)  frames 
respectively.  Details  of  coordinate  transformations  may  be  found  in  Chapter  2. 

£ej  is  the  elevation  angle  of  the  ejection  velocity  vector;  from  the  Z-model  s^j 
=  56.3°.  )Sej  is  the  azimuth  of  the  ejection  velocity  vector  and  is  measured  from  the 
north,  positive  clockwise  as  viewed  from  above  the  site  such  that  0  <  <  2tt,  and 

varies  from  (^  \  p  E[0:  nf6\  2ti)  for  a  total  of  13  possible  azimuths  to  | G  [0: 
11/36:  27t)  for  a  total  of  71  possible  azimuths  (specific  details  are  shown  in  figure 
captions  below). 

3.3.3.  Mass- Velocity  Distribution 

Eor  gravity-dominated  cratering,  volume  ejected  faster  than  a  given  velocity  is 
given  by  Equation  2.16,  where  Vgj  is  the  ejection  velocity  and  R  is  the  final  radius  of 
the  crater  [Holsapple,  1993;  Morrison  et  al,  2009].  Eor  the  coefficients  in  Equation 
2.16,  a  sand-like  surface  is  well  represented  [Melosh,  1989]  by  v  =  1.2,  and  mass- 
velocity  distributions  from  experimental  impacts  into  granular  targets  [Housen  and 
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Holsapple,  2011]  find  Cgy  =  0.25,  a  value  that  correlates  well  to  other  published 
values  between  0.25  and  0.36  [Andrews,  1975;  Stdjfler  et  al,  1975;  Cintala  et  al, 
1999].  We  adopt  Cgj  =  0.3  and  v  =  1.2.  Using  Equation  (2.16)  and  a  Phobos  density 
of  1.9  g/cc  [Avanesov  et  al.,  1989;  Rosenblatt  et  al.,  2008],  the  mass  ejected  within 
each  velocity  bin  is  calculated.  For  the  Stickney  impact,  of  the  1.5xl0'^  kg  ejected 
faster  than  the  Phobos  escape  velocity,  only  1.2xl0'^  kg  (8%)  is  ejected  faster  than 
100  ms  '.  Similarly,  with  2.72xl0'^  kg  ejected  faster  than  30  ms  ',  80%  of  ejecta  mass 
with  sufficient  velocity  to  escape  Phobos  leaves  between  11.3-30  ms  '.  Numbers  are 
nearly  identical  for  craters  down  to  1  km  diameter  formed  in  the  gravity  regime. 

3.4.  Results 

3.4.1.  Primordial  versus  modern  orbit 

We  find  that  a  majority  of  the  mass  ejected  from  Phobos  at  low  velocity  (<100 

ms’')  reaccretes;  the  rest  impacts  Mars  or  leaves  the  system.  A  minor  fraction  (<1%) 

impacts  Deimos  at  randomized  locations.  This  work  first  characterizes  Stickney,  the 

largest  crater  on  Phobos,  which  either  formed  geologically  recently  or  when  Phobos 

was  in  a  more  distant  orbit  close  to  the  synchronous  \mQ[Bums,  1978;  Cazenave  et 

al.,  1980;  Yoder,  1982].  Figure  3.1(a)  shows  the  reimpact  map  for  Stickney  had  it 

occurred  in  Phobos’  present  orbit,  with  a  mean  reaccretion  time  of  ~22  hours  where  a 

~  2.77  i?Mar5.  Catenae  from  Stickney  are  clearly  distinguishable.  Figure  3.1(b)  shows 

the  reimpact  map  for  a  primordial  a  =  6i?mars  [Yoder,  1982].  Here,  with  a  mean 

reaccretion  time  of  ~32  years  and  max  of  88  years,  individual  reaccreted  particles 
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have  lost  their  geometrical  association  due  to  the  much  longer  flight  time.  A  study  of 
sesquinaries  originating  from  Deimos  finds  that  ejecta  particles  reaccrete  over  -500 
years  and  do  not  form  crater  chains  \Nayak  et  al,  2016b].  Therefore,  for  Phobos’ 
modern  orbit  inside  the  Roche  limit,  reimpacts  are  not  only  more  frequent  than  in  a 
primordial  orbit,  but  occur  over  shorter  timescales  (days)  and  are  recognizably 
coherent,  forming  chains. 

3.4.2.  Nature  of  reimpacing  ejecta 


The  pattern  of  reimpacts  as  pictured  against  increasing  ejection  velocity 
(Figure  3.2)  shows  that  the  catenae  originate  from  very  low-velocity  particles  just 
above  the  escape  velocity  of  Phobos.  At  ejection  velocities  of  <  25  ms  ',  nearly  all 
reimpacts  cluster  in  catenae-like  patterns;  above  30  ms  ',  a  few  times  Vesc,  reimpacts 
are  less  correlated  to  catenae-like  structures.  Studying  the  reimpact  velocities  (Figure 
3.3),  we  find  that  they  are  correspondingly  higher  than  Vesc  [Davis  et  al,  1981; 
Thomas,  1998]  but  sufficiently  low  that  craters  produced  are  expected  to 
morphologically  resemble  secondary  craters.  In  addition,  though  the  mean  reaccretion 
time  for  all  Phobos-impacting  particles  is  -22  hours,  the  large  majority  of  catena¬ 
forming  impacts  occur  in  -6-14  hours. 

As  such,  these  are  a  special  case  of  sesquinary  impacts  (dosquinaries). 
Although  the  impactors  temporarily  escape  the  gravitational  influence  of  Phobos, 
they  exhibit  behavior  that  may  be  best  described  as  intermediate  between  a  secondary 
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and  traditional  sesquinary  impactor.  Since  80%  of  the  total  mass  that  escapes  Phobos 
during  an  impact  is  ejected  at  speeds  of  <30  ms  '  (Section  3.3),  this  mass  influx  is  a 
mechanism  capable  of  significant  impact  to  the  surface  geology;  for  detailed 
characterization  of  catena  production  this  work  now  focuses  on  10-30  ms  ' 
streamlines. 

Next,  we  attempt  to  establish  if  the  catenae  are  a  frequent  phenomenon  by 
searching  for  a  lower  limit  of  impact  crater  diameter  that  results  in  the  creation  of 
catenae-like  structures.  We  generate  <30  ms  '  reimpact  maps  for  ejecta  from  primary 
craters  of  1-km,  3-km  and  5-km  in  diameter  (Stickney:  ~10  km  diameter).  As  seen  in 
Figure  3.7,  catenae  are  predicted  to  form  as  a  byproduct  of  sesquinaries  from  primary 
craters  at  least  as  small  as  1-km  in  diameter,  and  likely  smaller.  The  saturated  (at  least 
down  to  300  m  craters  [Veverka  and  Duxbury,  1977])  surface  of  Phobos  suggests  that 
the  creation  of  low-velocity,  clustered  linear  impact  structures  from  sesquinary  ejecta 
is  a  relatively  frequent  process;  catena  formation  is  approximately  as  frequent  as 
gravity-regime  crater  formation,  where  most  of  the  escaping  ejecta  mass  is  just  barely 
escaping  (<30  ms  ').  This  suggests  that  linear  chains  of  low- velocity  impact  structures 
are  a  relatively  frequent  process  on  modern  Phobos,  each  correlated  with  a  source 
crater. 

It  is  interesting  to  note  that  a  kilometer-sized  crater  in  the  strength  regime 
would  create  faster  ejecta,  consequently  producing  less  reaccreting  ejecta  [Asphaug 
and  Melosh,  1993]  and  few  (if  any)  crater  chains.  This  analysis  therefore  allows  one 
to  probe  the  dynamics  of  crater  ejecta  in  a  way  not  done  before,  for  instance  proving 
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that  catena-producing  craters  formed  in  the  gravity  regime  (in  deep  regolith)  and  are 
likely  not  much  older  than  Phobos’  current  orbit  (<  50  Ma).  Figure  3.4  illustrates  how 
orbital  ejecta  lingering  in  the  vicinity  of  the  Phobos  orbit  can  be  swallowed  up  in 
hemispherical  patterns  that  lead  to  chain  or  catena-like  reaccretion  patterns  (e.g. 
Figure  3.1).  Tracing  the  precise  orbital  history  of  multiple  particles  shows  that  ejected 
particles  transition  from  the  gravitational  influence  of  Phobos  to  Mars  for  a  period  of 
between  1-5  orbits  before  subsequent  impact,  reaccreting  before  they  have  time  to 
disperse.  The  mean  reaccretion  time  of  ~22  hours  is  ~3  times  the  orbital  period  of 
Phobos;  for  <100  ms'^  impact  velocity  events,  the  mean  reaccretion  time  is  -6-14 
hours,  or  -1-2  times  the  orbital  period  of  Phobos. 

The  higher  the  ejection  velocity  (Figure  3.3),  the  higher  the  corresponding 
scatter  is  in  the  locations  of  reaccretion  (Figure  3.2).  This  correlates  to  longer  flight 
times  and  greater  interactions  with  Mars-dominated  gravity,  as  opposed  to  the  low- 
velocity  impactors,  which  escape  Phobos-centered  gravity  but  do  not  stray  far  from 
the  orbit  of  the  satellite  before  reimpact.  Low-velocity  impactors  creating  the  catena 
appear  to  be  dosquinary  in  nature,  i.e.,  intermediate  between  typical  secondary  and 
sesquinary  impactor  behavior. 

Given  a  low-velocity  ejection  bracket  (11-30  ms'\  Figure  3.5),  the  azimuth  of 
ejection  controls  the  length  and  coherence  of  the  catenae.  From  the  ejection  azimuths 
for  reaccreting  particles,  it  can  be  seen  that  different  azimuths  result  in  different 
“sections”  of  the  catena  being  formed  (shape  legend.  Figure  3.5).  The  formation  of  an 
axisymmetric  cone,  with  ejecta  along  every  outbound  azimuth,  is  of  course  an 


49 


idealized  case.  In  reality,  depending  on  the  geometry  of  the  primary  crater  formation 
(and  resulting  streamlines),  catenae  may  range  from  highly  focused  to  “patchy”.  If  no 
particle  is  ejected  between  the  azimuths  of  (say)  30°-45°,  there  will  be  a 
corresponding  “gap”  in  the  catena.  The  location  of  this  gap  can  be  inferred  from  the 
respective  shape  legends.  Inversely,  given  a  suspected  catena  on  Phobos,  this  makes  it 
possible  to  infer  the  properties  of  the  primary  impact. 

3.4.3.  Grid  Search  Varying  the  Primary  Crater  Location 

Good  correlation  exists  between  models  of  tidal  stresses  induced  by  the  orbital 
decay  of  Phobos  toward  Mars  and  several  groove  families,  but  cannot  match  a 
number  of  grooves  across  the  surface  [Asphaug  et  al,  2015b].  Can  the  observed 
catenae  detailed  in  this  work  account  for  these  misfit  grooves?  Even  among  this 
subset  of  grooves,  a  wide  variety  of  orientations  are  seen  [Hurford  and  Asphaug, 
2015].  To  determine  if  sesquinary  catenae  can  be  responsible  for  all  these  varied 
orientations,  we  investigate  the  effect  of  ejecta  originating  from  sources  other  than 
Stickney. 

An  equidistant  longitudinal  and  latitudinal  grid  is  now  created,  ranging  from 
90°  N  to  90°  S  in  30°  increments  and  180°  W  to  180°  E  in  90°  increments,  for  a  total 
of  35  grid  points.  This  span  accounts  for  variations  between  the  leading  and  trailing 
apex  of  Phobos,  as  well  as  the  sub-  and  anti-Mars  points.  Since  it  has  already  been 
shown  that  the  size  of  the  origin  crater  has  no  effect  on  the  geometry  of  the  resulting 
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reaccretions,  the  velocity  distribution  from  a  3 -km  diameter  crater  is  chosen  as  the 
standard  distribution  and  “released”  from  each  grid  point. 

The  results  presented  in  this  chapter  show  that  the  longitude  of  the  primary 
crater  has  no  effect  on  reaccreting  catenae  for  impacts  in  both  the  southern 
hemisphere  (Figure  3.8)  and  northern  hemisphere  (Figure  3.9).  In  these  figures,  note 
that  we  have  restricted  ourselves  to  showing  results  for  ^  |  G  [rr:  11/36:  2it)  for 
clarity;  results  are  mirrored  for  |  G  [0:  7t/36:  it).  This  can  be  seen  by  comparing 
the  panels  of  Figure  3.1 1,  the  reimpact  map  for  the  Stickney  impact. 

While  longitude  of  the  primary  crater  has  little  effect,  the  latitudinal  location 
of  the  primary  impact  has  a  direct  relation  to  the  orientation  of  the  catenae.  Figure  3.5 
shows  catenae  resulting  from  primary  impacts  at  0°,  30°,  60°  and  85°  N  latitude.  An 
interesting  pattern  emerges:  polar  primary  impacts  create  horizontal  catenae,  while 
vertical  chains  result  from  equatorial  primary  impacts.  These  results  are  mirrored  in 
the  southern  hemisphere  (Figure  3.11).  The  range  of  possible  catenae  orientations, 
from  horizontal  to  vertical,  show  that  low-velocity  sesquinary  impactors  could  indeed 
match  the  orientation  of  those  grooves  that  do  not  fit  a  tidal  evolution  origin. 
Orientations  are  mirrored  across  primary  impacts  to  either  hemisphere,  which  can  be 
seen  by  comparing  Figure  3.5  (northern  hemisphere)  to  Figure  3.11  (southern 
hemisphere).  From  the  variation  between  orbital  geometries  at  the  time  of  the  primary 
impact,  the  location  of  the  catenae  can  shift  longitudinally,  depending  on  the  location 
of  the  primary  impact  along  Phobos’s  orbit  around  Mars. 
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3.4.4.  Tracing  catenae  back  to  a  source  crater 


Based  on  the  relationships  of  reaccreting  catenae  to  a  primary  crater  location, 
we  can  now  match  a  catena  to  its  source  crater,  not  previously  done  for  any  planetary 
body.  In  doing  so  we  can  constrain  the  ejecta  velocity  field  and  provide  knowledge  of 
the  material  properties  in  the  region  of  the  source  crater.  For  our  case  study,  we 
choose  one  of  the  more  prominent  Phobos  catenae,  shown  in  Figure  3.6(a). 
Previously  studied  in  the  literature  {Veverka  and  Duxbury,  1977]  and  hemispheric  in 
extent,  it  is  not  obviously  correlated  with  Stickney  or  any  other  crater,  is 
morphologically  similar  to  cratering  expected  for  ejecta  colliding  at  just  above  ~Vesc 
and  crosscuts  the  predicted  stress  field  for  tidal  grooves  [Asphaug  et  al,  2015b; 
Hurford  and  Asphaug,  2015].  This  makes  it  a  suitable  test  of  the  hypothesis  that 
sesquinary  catenae  can  match  those  grooves  that  do  not  fit  tidal  models. 

In  comparing  to  Figure  3.5(d),  the  highlighted  catena  is  similar  to  reimpacts 
predicted  for  a  near-polar  primary  source  crater,  narrowing  the  grid  search  to  above 
60°  N.  Sesquinary  ejecta  emanating  from  an  impact  at  the  2.6-km  diameter  crater 
Grildrig  (81°  N,  196°  W)  is  modeled,  and  a  very  close  match  to  the  observed  catena  is 
found  (Figure  3.6).  This  test  case  shows  that  reimpacting  slow  co-orbital  ejecta  can 
explain  previously  mysterious  features  not  well  explained  by  any  previously  proposed 
mechanism.  These  ejecta  would  be  proximal  to  the  crater  on  a  'normal-gravity'  body 
like  the  Moon,  but  on  Phobos,  they  get  pulled  apart  into  strands  before  re-impact  a 
matter  of  orbits  later. 
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3.4.5.  Changing  the  Z-Model  Ejection  Angle 


Equation  2.6  fixes  the  ejection  angle  for  the  velocity  distribution:  in  this  final 
subsection,  it  is  investigated  whether  this  affects  the  outcome  of  the  simulations 
presented  above.  Figure  3.12.  shows  the  reaccretion  map  for  a  simulation  in  which 
the  ejection  angle  of  the  outbound  ejecta  is  allowed  to  stochastically  vary  between 
=  45  °  and  £ej  =  65  ° ,  i.e.,  for  Z-model  numbers  between  3  <  z  <  4.14 .  This  is 
overplotted  on  top  of  results  for  a  simulation  which  uses  our  standard  values  of  £ej  = 
56.3°  and  Z  =  3.5.  The  similarity  between  the  results  shows  that  changing  ejection 
angle  or  Z  number  does  not  change  the  qualitative  results  presented;  catena-like 
formations  are  still  expected.  However,  the  ejection  angle  may  represent  a  constraint 
on  the  crater  ejection  mechanism;  varying  the  Z-number  may  further  improve  the  fit 
of  our  modeled  catena  to,  for  example,  ejecta  from  the  Grildrig  crater. 

3.5.  Discussion 

The  results  presented  in  this  chapter  support  Grildrig’ s  formation  in  the 
gravity  regime  and  in  deep  regolith.  Grildrig  does  show  raised  rims  corresponding  to 
gravity-regime  production,  is  reasonably  fresh,  and  of  suitable  diameter  to  produce 
the  observed  catena.  When  loose  material  is  ejected  at  velocities  just  exceeding 
escape  velocity,  self-gravity  of  the  material  becomes  a  factor  and  leads  to  clump 
formation,  as  seen  from  asteroid  family  and  binary  formation  simulations  [Durda  et 
al,  2007;  Walsh  et  al,  2008].  This  is  how  an  impact  to  deep  regolith  can  release  large 
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ejecta  fragments,  which  would  then  become  sesquinary  impactors.  This  is  also  why 
the  catenae  do  not  have  cleanly-discriminated  crater  forms,  but  rather  clumpy 
streamers  of  interconnected  features. 

The  rim  of  Grildrig  shows  craters  similar  in  size  to  the  catena  in  its  vicinity;  it 
is  likely  that  reimpacting  ejecta  may  have  formed  some  of  these.  Better  image  data 
from  a  dedicated  Phobos  mission  can  enable  studies  of  chronology  in  the  relative 
sense,  i.e.,  evaluate  the  hypothesis  that  craters  on  the  rim  of  Grildrig  are  coeval  with 
catena  craters  associated  with  sesquinary  ejecta  from  Grildrig.  Certainly,  if  sesquinary 
ejecta  is  indeed  the  source  of  the  indicated  crater  chain,  there  is  little  one  can  say 
about  crater  ages  on  Phobos  from  the  consideration  of  sub-km  scale  crater  densities; 
age  estimates  from  crater  counting  were  previously  used  [Thomas  et  al,  1978]  as  an 
argument  against  a  tidal  origin  [Soter  and  Harris,  1977]  for  grooves. 

In  conclusion,  catenae  on  Phobos  created  by  low-velocity  sesquinary 
impactors  are  persistent  across  the  range  of  longitudinal,  latitudinal,  orbital  and 
conjunctive  variations,  with  distinctive  resulting  geometries.  This  in  turn  implies  that 
catena  formation  is  approximately  as  frequent  as  gravity-regime  crater  formation, 
where  most  of  the  escaping  ejecta  mass  is  just  barely  escaping  (<30  ms'^).  For  such 
events,  the  fate  according  to  our  model  is  for  material  to  impact  in  a  linear  chain. 
Based  on  Figure  3.3,  we  expect  the  resulting  crater  morphology  to  be  similar  to 
secondary  craters  in  nature,  consistent  with  secondary  crater  chains  and  catenae  noted 
on  the  Moon.  The  direct  association  of  sesquinary  catenae  with  source  craters  is  an 
important  and  new  kind  of  planetary  data  set,  leveraging  a  1:1  correspondence  to 
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constrain  the  dynamical  and  geologic  history  of  Phobos  in  a  novel  way.  Without 
applying  any  morphological  criterion,  we  can  say  if  a  catena-like  feature  on  Phobos  is 
a  good  candidate  for  a  sesquinary  origin  by  asking  if  it  follows  the  distinct  trending 
direction  that  it  must  follow  if  it  were  sesquinary  (Figure  3.5).  Based  solely  on  an 
inspection  of  Figure  3.5,  initial  guesses  may  be  made  at  the  latitude  of  that  feature’s 
possible  source  crater;  a  grid  search  similar  to  that  performed  for  Grildrig  would  then 
narrow  possibilities  down  to  one  source  crater. 

Significant  craters  without  corresponding  catenae  might  need  to  have  formed 
when  Phobos  was  in  a  more  distant  orbit,  or  in  strength-controlled  impacts  with  faster 
ejecta.  Conversely,  catenae  without  corresponding  source  craters  may  have  formed 
when  Phobos  was  in  a  different  orbital  or  tidal  locking  geometry  with  respect  to 
Mars.  The  lack  of  correlation  of  major  catenae  with  our  predictions  for  Stickney 
(Figure  3.1)  suggest  that  Stickney  formed  when  Phobos  was  more  distant  from  Mars 
[Schmedemann  et  al,  2014].  For  example,  long  flight  times  for  particles  reaccreting 
to  Deimos  result  in  many  reimpacts  but  no  catenae  on  that  body  [Nayak  et  al,  2016b]. 
This  does  not  preclude  the  possibility  of  impact  craters  on  Phobos  that  are 
uncorrelated  Stickney  sesquinaries. 

Discriminating  sesquinary  catenae  from  tidal  stress-induced  fissures  [Asphaug 
et  al,  2015b],  or  catenae-like  features  formed  from  regolith  draining  into  fissures 
[Horstman  and  Melosh,  1989],  is  the  next  step  in  Phobos  surface  science.  For 
example,  sesquinary  catenae  may  exhibit  raised  rims  typical  of  impact  craters, 
whereas  catenae-like  features  formed  by  regolith  draining  would  not.  Similarly,  the 
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walls  of  drainage  pits  would  stand  at  the  angle  of  repose,  whereas  sesquinary  craters 
might  be  less  steep.  However,  given  the  generally  heterogeneous  quality  and 
resolution  of  Phobos  images,  and  the  lack  of  systematic  mapping  products  applied  to 
the  Mars  Express  data  sets  [Witasse  et  al,  2014],  discriminating  catenae  as 
sesquinary  or  otherwise  solely  on  the  basis  of  imaging  is  an  expansive  task,  rendered 
further  difficult  by  the  fact  that  there  may  be  multiple  mechanisms  for  linear  feature 
formation  at  play.  There  is  a  need  for  published  topographic  profiles  of  Phobos  that 
reliably  measures  the  slopes  of  Phobos  features  relative  to  the  satellite’s  effective 
gravity  (which  can  vary  measurably  across  the  surface  of  the  moon),  from  which  the 
distinct  origin  of  a  particular  catena  might  be  distinguished.  The  creation  of  these 
products  is  ongoing. 

Using  this,  future  work  will  start  with  systematic  and  objective  geomorphic 
measurements  of  the  directions,  slopes,  depths,  and  the  sizes  and  spacings  of  linear 
pits  and  craters,  in  order  to  classify  linear  features  on  Phobos  objectively  (e.g.  into 
secondary  ejecta-like  streamers  versus  fissures,  etc.),  analogous  to  [Morrison  et  al, 
2009]  but  applied  to  modern  Phobos  data  [Witasse  et  al,  2014].  These  classified 
features  will  then  be  compared  to  model  predictions,  especially  sesquinary  catena 
predictions  from  fresh-looking  source  craters  in  present  and  geologically  recent 
orbits,  groove  predictions  from  the  tidal  model  [Asphaug  et  al,  2015b]  and  any  other 
viable  mechanisms.  This  study  also  has  potential  implications  for  other  closely  bound 
satellite  systems,  especially  Pluto  and  Charon,  whose  tidal  locking  and  close 
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proximity  leads  to  the  possibility  of  similar  short-timescale  reaccretions  collecting 


preferentially  onto  facing  hemispheres. 


3.6.  Figures 
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Figure  3.1.  Reaccretion  map  for  the  Stickney  impact  in  Phobos  modern  and  primordial  orbits, 
(a)  Ejecta  from  Phobos  in  its  modern  orbit,  a  -  2.77Rinars;  multiple  linear  strings  of  reimpacts  are 
noted  in  both  hemispheres  (brown  and  blue  ellipses),  (b)  Ejecta  from  Phobos  in  its  primordial 
orbit  near  the  Mars  synchronous  line,  a  -  6Rmars;  similar  reaccretion  patterns  are  not  seen.  All 
ejection  velocities  range  from  11-100  ms'^,  azimuth  p  6  [0:  7r/6:  In). 
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Figure  3.2.  The  evolution  of  catenae  with  increasing  velocity  of  ejection  from  Phobos.  Maps  show 
reimpact  locations  for  sesquinary  impactors  with  ejection  velocity  greater  than  Phobos  escape 
velocity  and  lesser  than:  (a)  20  ms  \  (b)  25  ms'^,  (c)  30  ms'^  and  (d)  100  ms'^.  Catenae-like 
formations  are  evident  from  particles  ejected  at  lower  velocities.  As  ejection  velocity  increases 
past  30  ms'^,  reimpacts  become  less  correlated  to  catenae-like  structures.  The  source  of  the 
catenae  is  therefore  very  low-velocity  particles.  The  primary  impact  location  is  Stickney  and 
azimuth  p  6  [0:  Tt/6:  In). 
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Figure  3.3.  Impact  vs.  ejection  velocity  for  sesquinary  impactors.  Ejecta  corresponds  to 
simulations  in  all  panels  of  Figure  3.2.  The  majority  of  impacts  occur  at  a  velocity  comparable  to 
their  ejection  velocity,  i.e.,  between  11-30  ms'^.  However  a  small  fraction  experience  acceleration 
due  to  greater  interactions  with  Mars’  gravity  (inset);  these  have  a  correspondingly  longer 
accretion  time.  Based  on  these  results,  sesquinary  impact  craters  are  largely  expected  to  have  an 
appearance  similar  to  secondary  craters  noted  in  image  surveys  of  Phobos. 


Figure  3.4.  Visualization  of  the  orbital  history  of  randomly  selected  particles.  Shown  are  the 
trajectories  of  randomly  selected  projectiles  ejected  at  11-30  ms'^,  which  impact  over  a  relatively 
short  time  in  a  grouped  fashion.  As  Phobos  moves  along  its  orbit  (toward  the  top  of  the  figure, 
dashed  pink  line),  ejected  particles  in  the  vicinity  of  Phobos’  orbit  (crimson,  white,  green,  red. 
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yellow  lines)  are  reaccreted  to  the  satellite.  Dashed  circles  represent  planetocentric  impact 
points,  illustrating  how  hemispherical  catena  (see  Figure  3.5)  may  be  formed. 
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Figure  3.5.  Catenae  orientations  change  from  vertical  to  horizontal  depending  on  the  latitude  of 
the  primary  impact.  Maps  show  resultant  catenae  from  primary  impacts  on  Phobos  at  the  prime 
meridian  and  (a)  0°  N;  (b)  30°  N;  (c)  60°  N;  (d)  85°  N.  Resulting  orientations  are  the  mirror 
inverse  of  those  from  southern  hemisphere  impacts  (compare  to  Figure  3.11).  Variations  with 
orbital  geometry  at  primary  impact  can  be  seen,  namely,  when  Phobos  is  at  Mars  periapsis  (red), 
apoapsis  (blue),  halfway  between  periapsis  and  apoapsis  along  the  ascending  (purple)  and 
descending  node  (green).  Ejection  velocities  are  11-30  ms'^;  source  crater  is  of  3-km  diameter. 
Shapes  denote  ejection  azimuth  (legend);  p  6  [n:  Tr/36:  2n)  are  shown  (compare  to  Figure  3.10). 
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Figure  3.6.  Tracing  an  observed  catena  back  to  its  source  crater,  (a)  Spacecraft  image  of  Phobos 
(Photo  credit:  ESA/Mars  Express)  showing  the  observed  catena  of  interest  (red  arrows);  (b) 
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Reimpact  map  for  a  primary  impact  at  Grildrig,  azimuth  p  6  [0:  7r/36:  2ii)  rendered  in  3D. 
Relative  sizes  and  orientations  between  (a)  and  (b)  are  similar  and  may  be  correlated  from 
labeled  craters;  Flimnap  and  Skyresh  craters  are  in  shadow  in  (a).  Features  underlined  in  yellow 
in  (b)  are  on  the  opposite  hemisphere.  From  the  correlation,  the  highlighted  catena  likely 
originates  from  sesquinary  ejecta  from  Grildrig. 
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Figure  3.7.  Impact  of  primary  crater  size  on  resultant  catenae.  Maps  show  resultant  catenae 
from  velocity  distributions  centered  at  the  location  of  Stickney  but  with  Z-model  crater  sizes  of 
diameter  (a)  5  km;  (b)  3  km  and  (c)  1  km.  Reimpacting  particles  have  ejection  velocities  of  11-30 
ms'^;  no  change  is  evident  between  varying  crater  sizes  in  this  velocity  range.  Craters  as  small  as 
1  km  in  diameter  appear  capable  of  creating  catenae-like  structures.  Ejection  azimuth  p  is 
henceforth  restricted  (compare  to  Figure  3.10)  to  p  E  [tt:  7r/36:  2ir);  shapes  denote  ejection 
azimuth  of  reaccreted  particles  (legend,  lower  left).  Changes  with  orbital  geometry 
configurations  at  the  time  of  the  primary  impact  are  shown,  namely,  when  Phobos  is  at  Mars 
periapsis  (red),  apoapsis  (blue),  halfway  between  periapsis  and  apoapsis  along  the  ascending 
node  (purple)  and  the  descending  node  (green). 


61 


-135  -90  -45  0  45  90  135  180 

Longitude  (deg) 


Figure  3.8.  Impact  of  primary  crater  longitude  (southern  hemisphere)  on  resultant  catenae. 
Maps  show  resultant  catenae  for  a  primary  impact  at  latitude  30°  S;  longitude  ranges  between 
(a)  90°  W;  (b)  prime  meridian;  (c)  180°  E.  Large  changes  in  primary  impact  crater  longitude 
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have  no  effect  on  the  orientation  of  the  resulting  catenae  and  are  degenerate  with  orbital  phasing 
of  Phobos  around  Mars.  Reimpacting  particles  have  ejection  velocities  of  11-30  ms'^  and  p  e[Ti: 
71/36:  271).  Colors  and  shapes  denoting  the  ejection  azimuth  are  as  in  Figure  3.7. 
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Figure  3.9.  Impact  of  primary  crater  longitude  (northern  hemisphere)  on  resultant  catenae. 
Maps  show  resultant  catenae  for  a  primary  impact  at  latitude  60°  N;  longitude  ranges  between 
(a)  90°  W;  (b)  prime  meridian;  (c)  180°  E.  Large  changes  in  primary  impact  crater  longitude 
have  no  effect  on  the  orientation  of  the  resulting  catenae  and  are  degenerate  with  orbital  phasing 
of  Phobos  around  Mars.  Reimpacting  particles  have  ejection  velocities  of  11-30  ms'^  and  p  E[ti: 
71/36:  271).  Colors  and  shapes  denoting  the  ejection  azimuth  are  as  in  Figure  3.7. 
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Figure  3.10.  Comparing  resultant  catenae  with  azimuths  of  ejecta  release.  Maps  show  resultant 
catenae  from  the  Stickney  impact  (centered  at  1°  N,  46°  W)  from  particles  with  ejection  velocities 
from  11-30  ms'^.  Colors  are  as  in  Figure  3.7.  (Top)  Azimuths  of  ejected  particles  are  p  \  p  e[0: 
77/36:  277);  (bottom)  compares  catenae  for  azimuths  restricted  to  P  \  P  E[ti:  77/36:  277).  Shapes 
denoting  ejection  azimuth  are  in  legends  (bottom  left)  of  both  panels.  Mirroring  the  catenae  in 
the  bottom  panel  yields  the  full  picture;  for  clarity  in  comparing  geometry  configurations  (GCs) 
we  show  the  subset  of  azimuths  in  the  bottom  subfigure. 
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Figure  3.11.  Impact  of  primary  crater  latitude  (southern  hemisphere)  on  resultant  catenae.  Maps 
show  resultant  catenae  from  primary  impacts  on  Phobos  at  the  prime  meridian  and  (a)  0°  S;  (b) 
30°  S;  (c)  60°  S;  (d)  85°  S.  Catenae  orientation  can  change  from  near- vertical  to  horizontal, 
depending  on  the  latitude  of  the  primary  impact,  and  mirror  northern  hemisphere  impacts 
(compare  to  Fig  5).  Reimpacting  particles  have  ejection  velocities  of  11-30  ms'^,  p  6  [n:  n/ 
36:  271);  colors  and  azimuth  legend  as  in  Figure  3.7. 
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Figure  3.12.  Variation  of  resultant  catenae  with  Z-number  (ejection  angle).  Reaccretion  map  for 
ejecta  from  the  Grildrig  crater,  at  Phobos  periapsis  around  Mars,  where  the  ejection  angle  is 
held  fixed  at  56.3°  (blue  catena)  and  allowed  to  vary  stochastically  between  45°  and  65°  (purple 
catena).  Highly  similar  reaccretion  patterns  between  the  two  methods  show  that  the  catena 
formation  mechanism  detailed  here  does  not  depend  on  Z-number  or  ejection  angle  variations. 
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Chapter  4 


Magnetics:  A  Case  Study  of  the  Lunar 
South  Pole-Aitken  Basin 

This  chapter  is  a  modified  reprint  ofM.  Nayak,  D.  Hemingway  and  I.  Garrick-Bethell 
(2016),  Magnetization  in  the  South  Pole-Aitken  Basin:  Implications  for  the  lunar 
dynamo  and  true  polar  wander,  accepted,  Icarus. 

4.1  Abstract 

A  number  of  magnetic  anomalies  are  present  along  the  northern  edge  of  the 

lunar  South  Pole-Aitken  (SPA)  basin.  A  variety  of  hypotheses  for  their  formation 

have  been  proposed,  but  an  in-depth  study  of  their  properties  has  not  been  performed. 

Here  we  use  two  different  methods  to  invert  for  their  source  body  characteristics:  one 

that  completely  searches  a  small  parameter  space  of  fewer  than  ten  uniform- strength 

dipoles  per  anomaly,  and  another  that  uses  grids  of  hundreds  of  dipoles  with  variable 

magnetization  strengths.  Both  methods  assume  uniform  magnetization  directions  at 

each  anomaly  and  with  one  exception,  produce  nearly  the  same  results.  We  introduce 

new  Monte  Carlo  methods  to  quantify  errors  in  our  inversions  arising  from  Gaussian 

time-dependent  changes  in  the  external  field  and  the  uncertain  geometry  of  the  source 

bodies.  We  find  the  errors  from  uncertainty  in  source  body  geometry  are  almost 

always  higher.  We  also  find  a  diverse  set  of  magnetization  directions  around  SPA, 

which  we  combine  with  other  physical  arguments  to  conclude  that  the  source  bodies 
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were  likely  magnetized  in  a  dynamo  field.  Igneous  intrusions  are  a  reasonable 
explanation  [Purucker  et  al,  2012]  for  the  directional  variability,  since  they  could  be 
intruded  over  different  magnetic  epochs.  However  the  directional  variability  also 
implies  either  surprisingly  large  amounts  of  true  polar  wander  or  non-axially  aligned 
dynamo  fields.  We  also  explore  the  possibility  that  true  polar  wander  caused  by  the 
SPA  impact  could  allow  iron-rich  SPA  ejecta  to  record  a  diverse  set  of  magnetic  field 
directions.  Some  of  this  material  may  have  also  become  “sesquinary”  ejecta  and  re¬ 
impacted  across  the  Moon  on  10"^- 10^  year  timescales  to  capture  such  changes.  No 
completely  satisfactory  answer  emerges,  except  that  the  dipole-axis  of  the  lunar 
dynamo  may  have  been  variable  in  direction. 

4.2  Introduction 

The  South  Pole-Aitken  (SPA)  Basin  is  the  Moon’s  largest  and  oldest  well- 
defined  basin  [Garrick-Bethell  and  Zuber,  2009].  It  is  also  the  site  of  the  largest 
grouping  of  magnetic  anomalies  on  the  Moon  [Purucker,  2008].  Initially  discovered 
by  the  Apollo  15  subsatellite  [Coleman  et  al,  1972],  the  origin  of  these  features 
remains  unknown.  Across  the  last  20  years,  three  dominant  hypotheses  for  their 
formation  have  emerged.  The  first  is  that  ejecta  from  the  Imbrium  and  Serenitatis 
basins  has  accumulated  at  their  antipodes,  which  fall  within  the  SPA  basin  (Hood  and 
Huang,  1991;  Lin  et  al.,  1998;  black  crosses  in  Figure  1).  This  antipodal  ejecta  then 
attains  a  remanent  magnetization  from  either  a  dynamo  field,  or  a  field  related  to  the 
impact  event  [Hood  and  Artemieva,  2008].  The  second  hypothesis  is  that  the 
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anomalies  arise  from  subsurface  dikes  that  cooled  in  a  dynamo  field  [Purucker  et  al, 
2012],  The  third  is  that  iron-rich  material  derived  from  the  SPA  impactor  cooled  in  a 
dynamo  field  [Wieczorek  et  al,  2012],  Determining  which,  if  any  of  these  hypotheses 
are  true,  would  have  implications  for  understanding  the  nature  and  history  of  the  lunar 
dynamo. 

Despite  the  importance  of  these  anomalies  in  understanding  lunar  magnetism, 
no  detailed  studies  of  the  source  body  characteristics  have  been  performed.  Purucker 
et  al.  (2012)  modeled  one  of  the  anomalies  within  SPA  as  a  series  of  vertically 
magnetized  dikes,  and  estimated  their  magnetization  strength.  Global  studies  of  lunar 
magnetic  anomalies  have  neglected  the  SPA  region  [Arkani-Hamed  and  Boutin, 
2014;  Takahashi  et  al,  2014],  because  of  the  complicated  field  structure  in  the 
region.  Using  a  method  they  refer  to  as  surface  vector  mapping,  Tsunakawa  et  al. 
(2014)  calculated  the  statistics  of  the  declination  of  the  field  over  the  anomalies  at 
SPA.  Based  on  the  distribution  of  declinations,  they  found  evidence  that  the  source 
bodies  are  horizontally  elongated  in  the  east-west  direction.  While  they  also  estimated 
the  magnetization  direction  and  magnetic  paleopole  from  the  anomaly  centered  on  the 
Leibnitz  crater  within  SPA,  they  did  not  estimate  the  source  body  characteristics  of 
any  of  the  larger  magnetic  anomalies  that  characterize  the  region. 

Here,  we  will  show  that  many  individual  anomalies  can  be  identified  within 
SPA.  We  use  software  that  can  readily  generate  and  compare  magnetic  field  maps 
from  all  available  observations,  helping  us  avoid  spurious  or  poorly  defined  magnetic 
anomalies  that  may  have  complicated  other  efforts  to  study  the  region.  Comparisons 
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with  a  spherical  harmonic  model  of  the  magnetic  field  \Purucker  and  Nicholas,  2010] 
provide  a  further  test  of  consistency.  We  will  show  that  the  source  body 
magnetization  directions  offer  a  key  test  for  their  formation  hypotheses,  and  we  make 
a  substantial  effort  to  characterize  the  uncertainty  in  these  directions.  The  paper  is 
organized  as  follows:  Section  4.3  presents  our  methods,  including  uncertainty 
estimation.  Section  4.4  presents  our  results,  and  in  Section  4.5,  we  discuss  possible 
origins  of  the  observed  diversity  in  magnetization  directions.  We  consider  an  impact 
versus  dynamo  origin  for  these  anomalies,  secular  variation,  dynamos  misaligned 
with  the  lunar  spin  axis,  and  both  long  and  short-timescale  true  polar  wander.  Finally, 
Section  4.6  contains  our  conclusions. 

4.3  Methods 

4.3.1  Data  Sources 

We  use  magnetometer  data  from  two  independent  sources:  Lunar  Prospector 
(LP-MAG)  and  SELENE/Kaguya  (K-MAG).  The  Lunar  Prospector  fluxgate 
magnetometer  measured  the  vector  magnetic  field  at  up  to  18  Hz  and  transmitted  its 
measurements  at  a  reduced  resolution  of  9  Hz.  Level  IB  LP-MAG  data  are  obtained 
from  the  NASA  Planetary  Data  System  (ppi.pds.nasa.gov).  The  SELENE/Kaguya 
spacecraft  used  a  tri-axial  fluxgate  magnetometer  with  a  sampling  rate  up  to  32  Hz. 
K-MAG  magnetometer  data  are  obtained  from  the  SELENE  data  archives 
(12db.selene.darts.isas.jaxa.jp).  The  cadence  of  measurements  used  in  this  study  is  0.2 
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Hz  for  LP-MAG  and  0.25  Hz  for  K-MAG.  Topography  data  are  from  the  Lunar 
Orbiter  Laser  Altimeter  (LOLA)  [Smith  et  al,  2010]  (pds-geosciences.wustl.edu). 

To  best  capture  the  Moon’s  crustal  field,  all  data  used  for  analysis  were 
collected  in  either  the  lunar  wake  or  while  the  Moon  was  in  the  Earth’s  magnetotail 
(wake/tail),  avoiding  distortions  caused  by  the  solar  wind  noted  by  [Kurata  et  al, 
2005;  Halekas  et  al,  2008].  Tail  datasets  specifically  exclude  epochs  during  which 
plasma  sheet  disturbances  were  noted  [Halekas  et  al,  2012].  Consecutive  orbits  are 
~1°  in  longitude  apart.  At  0.2  Hz,  successive  magnetometer  measurements  are 
separated  by  ~8  km  in  the  latitudinal  direction.  All  data  used  in  this  study  are  from 
the  final  months  of  Lunar  Prospector  in  1999  and  Kaguya  in  2009,  when 
measurements  were  taken  at  observation  altitudes  below  50  km. 

4.3.2  Data  Processing  and  Anomaly  Identification 

Details  on  the  generation  of  magnetic  field  maps  used  in  this  paper  (e.g. 
Figure  4.8)  are  after  Hemingway  and  Garrick-Bethell  (2012),  summarized  here  for 
completeness.  After  subtracting  the  background  field  (taken  to  be  the  mean  field 
across  each  orbit  segment  spanning  the  region  of  interest),  the  remaining  fields  are 
assumed  to  be  due  to  crustal  sources.  Data  from  consecutive  orbits  are  combined  and 
fit  to  square  meshes  using  Delaunay  triangulation,  with  a  grid  spacing  of  0.25°  x 
0.25°  (7.6  km  x  7.6  km  equatorial).  This  spacing  is  finer  than  the  spacing  between 
observations  (e.g.,  ~8  km  latitude,  ~30  km  longitude  for  LP-MAG),  ensuring  no  loss 
of  signal  variation  during  grid  generation.  Orbit  segments  that  visually  appear 
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distorted  by  transient  magnetic  noise,  despite  these  efforts,  are  discarded.  Gaps  in 
spacecraft  coverage,  if  any,  can  be  seen  in  magnetic  field  maps,  as  the  spacecraft 
measurement  locations  are  shown  in  all  figures  (e.g.  Figure  4.8,  white  points). 

To  analyze  the  northern  SPA  region,  we  divide  the  strongest  magnetic 
anomalies  into  ten  study  areas  (Figure  4.1).  We  choose  anomalies  whose  Kaguya  and 
LP  magnetic  field  maps  are  consistent  with  maps  from  a  spherical  harmonic  model  of 
the  field  [Purucker  and  Nicholas,  2010],  which  can  be  seen  by  comparing  Figure  4.1b 
and  Ic.  Anomalies  that  were  not  consistent  with  this  model  or  showed  artifacts  of 
external  disturbances  were  not  used.  Areas  1  and  2  encompass  two  approximately 
linear  and  perpendicular  magnetic  features;  in  their  vicinity  is  Area  3,  also  in 
proximity  to  the  Van  de  Graaff  crater  [Dyal  et  al,  1974].  Area  1  is  very  similar  to  the 
area  modeled  by  Purucker  et  al.  (2012).  Together,  we  refer  to  Areas  1-3  as  the 
northwestern  cluster.  Area  9  is  to  the  southwest  of  these  areas  and  is  associated  with 
the  “swirl”  albedo  anomalies  at  Mare  Ingenii  [Blewett  et  al,  2011;  Kramer  et  al, 
2011].  Area  4  is  an  isolated  anomaly  to  the  southeast  of  areas  1-3.  Areas  5-8  are  east- 
west  trending  linear  features  found  to  the  east  of  these  study  areas,  and  together  we 
refer  to  them  as  the  eastern  stripes.  The  anomalies  at  areas  1  and  8  appear  connected 
in  some  maps  of  the  magnetic  field  (Purucker  et  al.  2012).  Finally,  we  also  define 
area  10,  interesting  because  it  falls  outside  the  basin’s  outer  topographic  rim 
[Garrick-Bethell  and  Zuber,  2009]. 

We  divide  all  data  in  this  paper  by  the  day  of  observation,  such  that  each 
dataset  for  a  given  anomaly  is  collected  at  approximately  the  same  altitude.  The 
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datasets  may  be  referenced  to  the  day  on  which  they  were  collected,  e.g.,  1999  day 
172  (LP-MAG;  read  as  99172)  or  2009  day  123  (K-MAG;  read  as  09123).  Figure  4.2 
shows  an  overview  of  the  total  magnetic  field  maps  created  for  all  observation  days 
used  for  all  study  areas;  Figure  4.23  -  Figure  4.57  show  details  of  magnetic  inversions 
for  all  these  areas  (see  Section  4.3.3  and  4.3.4). 

4.3.3  Inversion  Algorithm  1:  Defined  Dipoles,  Constant 
Magnetization 

To  assess  the  robustness  of  our  results,  we  use  two  different  algorithms  to 
invert  for  source  body  characteristics.  Both  are  regressions  to  find  the  least  squared 
error.  The  first  completely  searches  a  given  parameter  range  using  manually  placed 
dipoles  as  the  source  bodies.  The  regression  parameter  space  varies  burial  depth  in 
km,  magnetic  dipole  moment  in  Am^,  magnetic  dip  (inclination)  in  degrees 
downward  from  the  horizontal,  and  declination  in  degrees  clockwise  from  north.  All 
dipoles  in  a  given  study  area  are  constrained  to  be  at  the  same  depth,  moment  and 
direction.  Source  dipole  positions  (latitude  Vj  and  longitude  0s)  placed  manually 
by  the  user  at  locations  where  the  magnetic  field  strength  is  greatest  according  to  our 
maps.  We  refer  to  this  henceforth  as  the  Defined  Dipoles,  Constant  Magnetization 
(DD-CM)  algorithm.  As  we  will  show,  the  ability  to  completely  search  a  parameter 
space  at  relatively  fast  computation  speeds  gives  this  algorithm  some  advantages 
when  uncertainties  must  be  estimated. 
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The  depth  is  allowed  to  vary  between  1  km  and  99  km  in  steps  of  0.25  km, 
and  the  dipole  magnetic  moment  between  lO'^  Am^and  10'"^Am^in  steps  of  2.5x10^' 
Am  .  Depth  solutions  are  further  constrained  to  be  no  shallower  than  the  deepest 
negative  topography  in  a  given  study  area,  ensuring  that  all  solutions  lie  beneath  the 
lunar  surface;  depth  is  measured  against  a  reference  sphere  of  1737.4  km  [Smith  et 
al,  2010].  The  magnetic  inclination  is  allowed  to  vary  from  -90°  to  -1-90°  in  steps  of 
1°  and  declination  from  -180°  to  -1-180°  in  steps  of  1°. 

For  n  field  measurements  inside  a  given  study  area,  the  difference  between  the 
model  and  the  data  is  computed.  We  minimize  the  mean  of  the  root  mean  square 
(RMS)  total  error  of  the  east  (SBgast)^  north  (SEfi^rth)  radial  (SBradiai) 
component  errors,  5Stotai’ 

total)  RMSS  ~  ^~^i=iC^^east  T  ^^north  (  4.1  ) 

This  quantity  is  then  ranked  in  decreasing  order  of  total  error  to  find  the  best- 
fit  magnetic  characteristics. 

4.3.4  Inversion  Algorithm  2:  Gridded  Dipoles,  Variable 
Magnetization 

We  test  the  validity  of  our  inversion  results  using  an  alternative,  more 

objective  approach,  at  the  expense  of  exploring  a  much  larger  parameter  space. 

Instead  of  placing  dipoles  manually,  we  establish  a  0.25°  x  0.25°  grid  of  dipoles 

covering  the  entire  study  area  [Nicholas  et  ah,  2007;  Hemingway  and  Garrick- 

Bethell,  2012].  All  dipoles  are  again  constrained  to  be  at  the  same  depth,  inclination 

73 


and  declination.  Unlike  the  first  algorithm,  however,  the  magnetic  moment  is  allowed 
to  vary  among  the  dipoles.  This  creates  a  large  n+3  parameter  space,  where  n  is  the 
number  of  dipoles;  Table  4.1  lists  the  number  of  gridded  dipoles  for  each  study  area, 
which  ranges  from  320-1408.  The  solution  is  then  found  via  a  genetic  search 
algorithm  that  minimizes  the  RMS  error  (for  details,  see  Hemingway  and  Garrick- 
Bethell  (2012)).  The  algorithm  evolves  through  “generations”  to  progress  toward  a 
better  fit  (smaller  error)  to  the  data.  However,  the  evolutionary  nature  of  the 
algorithm  does  not  guarantee  optimality  of  the  solution.  We  refer  to  this  henceforth  as 
the  Gridded  Dipoles,  Variable  Magnetization  (GD-VM)  algorithm.  The  grid  of  best- 
fit  magnetic  moments  found  by  this  method,  ranging  across  several  orders  of 
magnitude  at  each  anomaly,  can  be  seen  in  Figure  4.58  -  Figure  4.67,  which  show 
representative  results  for  each  study  area. 

4.3.5  Uncertainty  Estimation  for  Magnetization  Directions 

Error  in  our  regressions  arises  from:  1 )  time-variable  contributions  to  the 
measured  field  due  to  non-crustal  fields  or  instrument  noise,  2)  the  ideal  assumption 
that  all  sources  are  uniformly  magnetized  in  the  same  direction,  and  3)  the  simplified 
geometry  of  our  source  models,  even  if  the  assumption  of  unidirectional 
magnetization  were  completely  true. 

4.3.5. 1  Time  variable  contributions 

One  method  to  account  for  the  first  source  of  error  is  to  report  the  set  of 
solutions  with  an  RMS  error  equal  to  or  less  than  the  uncertainty  in  the  magnetic  field 
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measurements  [Parker,  1991].  The  range  of  magnetization  directions  in  this  set  then 
defines  the  error  in  the  best-fit  direction.  The  drawback  of  this  method  is  that  the  use 
of  the  field’s  uncertainty  is  arbitrary.  For  example,  at  different  observation  altitudes, 
the  magnitude  of  the  RMS  error  for  the  best-fit  solution  will  vary,  such  that  the 
measurement  uncertainty  will  have  a  different  impact  on  the  size  of  the  set  of 
solutions  that  represent  the  error.  More  importantly,  as  we  show  below,  the  set  of 
allowable  solutions  suggested  by  Gaussian  measurement  noise  (or  external  field 
oscillations)  of  a  given  magnitude  is  often  small  compared  with  the  size  of  the  set  of 
solutions  contained  by  an  RMS  error  of  equal  magnitude.  That  is,  using  the 
uncertainty  in  the  field  measurement  to  define  the  allowable  solution  set  can 
overestimate  the  error  in  the  regression,  at  least  in  the  case  of  lunar  magnetic  field 
measurements. 

To  demonstrate  this  effect,  we  simulate  the  effect  of  Gaussian  measurement 
noise  in  each  of  the  45  regressions  we  perform  (Figure  4.2).  Before  doing  so,  we  first 
estimate  the  characteristic  noise  magnitude.  For  both  LP-MAG  and  K-MAG,  the 
instrument  noise  is  <  0.1  nT,  such  that  the  dominant  source  of  error  arises  from 
fluctuations  in  the  interplanetary  magnetic  field  (IMF)  rather  than  the  instrument.  To 
estimate  a  characteristic  amplitude  of  the  oscillations  in  the  IMF  over  the  timescale  of 
observation  at  one  magnetic  anomaly,  while  in  the  wake  or  magnetotail,  one  would 
ideally  examine  the  time-dependent  oscillations  just  prior  to  flying  over  the  anomaly 
of  interest.  However,  there  is  no  simple  way  to  separate  oscillations  arising  from 
small-scale  crustal  fields  and  time  dependent  fields.  As  a  proxy,  we  can  examine  the 
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field  variations  over  the  Mare  Imbrium  region  while  the  Moon  is  in  the  wake.  Since 
the  Imbrium  region  is  known  to  have  the  lowest  amount  of  crustal  magnetism 
anywhere  on  the  Moon  [Mitchell  et  al,  2008;  Tsunakawa  et  al,  2014],  observations 
there  are  the  most  likely  to  represent  only  the  time  dependent  oscillations.  We  find 
that  the  standard  deviations  at  Imbrium  are  0.07  nT  and  0.13  uT  for  two 
representative  days  (Figure  4.4).  However,  to  be  conservative,  we  use  an  order  of 
magnitude  higher  standard  deviation  of  1  nT  as  the  representative  value.  Takahashi  et 
al.  (2014)  make  a  similar  (1  nT)  assumption  about  the  magnitude  of  external  field 
contributions. 

Next,  to  simulate  the  effect  of  typical  time-dependent  fields  in  our  DD-CM 
regressions,  we  added  random  Gaussian  noise  with  a  standard  deviation  of  1  nT  to 
each  field  component  (east,  north,  radial)  for  all  45  datasets.  An  example  is  shown  in 
Figure  4.5.  Note  that  we  are  assuming  that  the  noise  at  each  measurement  point  is 
uncorrelated  with  the  previous  measurement  point.  Our  assumption  is  only  valid  if 
oscillations  in  the  IMF  occur  on  timescales  faster  than  the  time  between 
measurements  (5  seconds  for  LP  data,  and  4  seconds  for  Kaguya  data).  Future  work 
will  determine  the  dependence  of  our  error  estimates  on  a  variety  of  correlation 
timescales,  and  better  determine  the  actual  timescale  of  IMF  oscillations. 

For  each  of  the  45  datasets  (Figure  4.2),  we  generated  20  sets  of  noise-added 
data.  Assuming  the  best-fit  directions  are  Fisher  distributed  [Fisher,  1953],  their 
angular  standard  deviations  and  Fisher  distribution  precision  parameters  k,  were  then 
obtained  from  each  set  of  20  simulations.  The  mean  k  value  across  all  datasets  was 
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then  used  as  the  final  estimate  of  the  uncertainty  from  time-dependent  fields.  The 
ability  to  simulate  the  effect  of  noise  in  this  manner  represents  an  advantage  of  the 
DD-CM  method,  which  has  a  faster  computation  time  compared  to  the  GD-VM 
method  (see  Figure  4.5). 

We  find  that  the  angular  standard  deviation  in  best-fit  direction  associated 
with  1  uT  noise  is  always  smaller  than  the  set  of  directions  permitted  when  using  an 
error  threshold  of  1  nT.  For  example,  the  best-fit  RMS  error  for  Area  1,  09096 
(Figure  4.5)  is  2.4  nT,  which  is  already  higher  than  the  ~1  nT  uncertainty  in  the  field 
(this  is  true  for  all  anomalies  in  this  study).  The  set  of  solutions  that  are  allowed  by 
considering  1  nT  additional  RMS  error  (total  residual  error  3.4  nT)  produces  -30°  of 
uncertainty  in  the  best-fit  magnetization  direction  (Figure  4.6).  In  contrast.  Figure  4.6 
also  shows  that  adding  1  nT  Gaussian  noise  to  these  data  via  the  above  method  only 
produces  an  angular  standard  deviation  of  5°  in  the  best-fit  direction.  The  smaller  size 
of  the  effect  of  simulated  Gaussian  noise  is  true  in  all  of  our  regressions. 

4.3.5.2  Non-uniform  magnetization  directions 

The  second  contribution  to  error  is  the  non-uniform  magnetization  directions 
in  the  source  bodies.  However,  it  is  impossible  to  make  any  inferences  about  the 
magnetization  if  we  permit  the  infinite  number  of  possible  source  magnetizations 
with  mixed  directions.  Therefore  we  must  at  least  assume  that  the  source  is  uniformly 
magnetized.  We  attempt  to  mitigate  this  effect  by  selecting  small,  well-defined 


77 


anomalies  (Figure  4.1).  Of  course,  if  non-uniformity  dominates  the  source  of  error 
then  these  regression  results  are  less  meaningful. 

4.3.53  Source  body  geometry 

Finally,  we  address  what  is  likely  the  dominant  source  of  uncertainty:  the 
complex  geometry  of  the  magnetic  source  bodies  compared  to  our  simplified  models. 
We  do  this  in  two  ways.  The  first  is  to  use  compute  the  best-fit  magnetization 
separately  for  different  altitude  datasets  at  each  anomaly.  Data  from  different  altitudes 
display  a  diversity  of  magnetic  field  morphologies,  as  field  strength  decays  with 
altitude  (which  ranges  in  this  study  from  13-46  km).  Hence,  each  altitude  will  lead  to 
different  choices  of  dipole  locations;  for  example,  compare  Figure  4.8  to  Figure  4.23  - 
Figure  4.25.  Similarly,  the  geometries  of  the  best-fit  source  bodies  will  also  vary. 
This  variability  is  an  advantage,  as  it  allows  us  to  probe  the  sensitivity  of  our  results 
to  source  geometry,  with  the  spread  in  the  best-fit  magnetization  directions 
representing  the  measure  of  uncertainty.  This  is  analogous  to  the  practice  in 
paleomagnetism  of  sampling  a  single  rock  formation  multiple  times  at  different  sites 
[Irving,  1964].  Because  each  measured  direction  at  each  altitude  is  independently 
calculated  we  can  also  assign  a  95%  confidence  interval  to  the  mean  direction  [Butler, 
1998].  In  practice,  we  combine  the  results  from  both  the  DD-CM  and  GD-VM 
algorithms  to  calculate  the  final  mean  and  confidence  interval,  even  though  the 
direction  measurements  for  a  given  day  are  not  strictly  independent  across  the  two 
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algorithms.  We  also  include  one  best-fit  direction  obtained  from  the  merged  data  at 
all  altitudes,  using  the  GD-VM  algorithm. 

The  second  way  we  account  for  uncertainty  due  to  source  body  geometry  is  to 
explore  the  range  in  directions  that  are  returned  by  displacing  the  nominal  DD-CM 
algorithm  dipoles  for  all  45  datasets.  We  calculate  the  effect  of  placing  dipoles 
anywhere  randomly  on  a  0.5°-radius  circle  from  their  nominal  location  (Figure  4.7). 
The  value  of  0.5°  represents  the  approximate  error  in  longitudinal  uncertainty  in  the 
anomaly  peak  field  location  (orbits  are  spaced  by  approximately  1°  for  a  given 
constant  altitude  data  set).  For  each  data  set,  we  use  100  random  placements  of 
dipoles  at  any  location  on  the  defined  circles  and  calculate  the  100  best-fit 
magnetization  directions  and  k  values.  The  k  values  are  then  averaged  within  an  area 
to  obtain  the  uncertainty  from  source  geometry  for  that  area.  Again,  this  represents 
an  advantage  of  the  DD-CM  method  compared  to  the  GD-VM  method,  due  to  the 
ability  to  modify  model  dipole  placements  and  a  faster  computation  time.  We  will  see 
this  advantage  manifest  when  exploring  ambiguous  regression  results  for  area  2. 
Finally,  we  note  that  for  area  4  we  use  0.25°-radius  circles,  due  to  the  close  proximity 
of  the  model  dipoles  there. 

The  dispersions  in  best-fit  directions  obtained  from  our  Monte  Carlo 
simulations  for  time-dependent  fields  and  0.5°-radius  dipole  displacements  are  then 
combined  using  their  k  values,  averaged  together  for  each  area  [Irving,  1964]: 


total 


time -dependent 


source -geometry 


(4.2) 


An  angular  standard  deviation  ^53  is  then  estimated  by  (Irving,  1964): 
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Finally,  it  is  added  to  the  angular  95%  confidence  interval  obtained  from  the 
variable  altitude  results  within  a  given  area. 

The  following  is  a  summary  of  our  uncertainty  estimation  methods: 

1 .  Account  for  time-dependent  external  field  contributions  and 
instrument  noise:  Perform  Monte  Carlo  simulations  for  every  anomaly  and 
altitude  with  the  addition  of  1  nT  noise,  and  obtain  mean  k  for  the  best-fit 
directions  within  an  anomaly,  using  the  DD-CM  algorithm  (Figure  4.5). 

2.  Account  for  uncertainty  in  source  geometry:  Perform  Monte 
Carlo  simulations  for  every  anomaly  and  altitude  while  altering  the  dipole 
placement,  and  obtain  mean  k  for  the  best-fit  directions  within  an  anomaly 
using  the  DD-CM  algorithm  (Figure  4.7,  Table  4.2). 

3.  Combine  the  k  values  from  steps  1  and  2  and  obtain  the  angular 
standard  deviation  via  equations  2  and  3. 

4.  Account  for  uncertainty  in  source  geometry  by  using  variable 
altitude  data,  and  hence  variable  dipole  placement.  Obtain  a  95% 
confidence  interval  on  the  mean  direction  from  the  combined  best-fit 
directions  from  the  DD-CM  and  GD-VM  algorithms  (Table  4.3). 

5.  Add  the  angular  dispersion  from  step  3  to  the  angular  95% 
confidence  interval  from  step  4  to  obtain  the  total  angular  uncertainty 
estimate  in  the  best-fit  direction. 
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4.4  Results 


4.4.1  Regression  Results 

Minimum  error  magnetization  directions  are  found  for  each  Lunar  Prospector 
and  Kaguya  dataset,  using  both  the  DD-CM  (Section  4.3.3)  and  the  GD-VM  (Section 
4.3.4)  algorithms.  Table  4.1  compiles  the  values  of  these  best-fit  directions.  The  best- 
fit  models  are  shown  in  Figure  4.8  -  Figure  4.17  for  study  areas  1-10.  The  figures 
show  10  representative  examples,  one  for  each  study  area;  Figure  4.23  -  Figure  4.57 
show  all  results  for  every  dataset  listed  in  Table  4.1.  Angular  standard  deviations  and 
k  for  the  noise-added  and  dipole  displacement  simulations  are  shown  in  Table  4.2  for 
all  datasets.  RMS  error  maps  for  all  datasets  (similar  to  Figure  4.6)  are  compiled  in 
Figure  4.18  and  Figure  4.19;  these  illustrate  the  difference  between  using  an  arbitrary 
uncertainty  threshold  and  the  Monte  Carlo  methods  employed  in  this  work. 

Figure  4.20  compiles  the  values  of  the  best-fit  directions  from  the  DD-CM 
and  GD-VM  algorithms.  For  the  DD-CM  algorithm,  we  show  one  standard  deviation 
of  dispersion  from  1  nT  noise  (Figure  4.20a)  and  one  standard  deviation  of  dispersion 
from  displacing  the  model  dipoles  by  0.5°  (Figure  4.20b)  (see  Methods).  Overall,  we 
find  the  error  from  the  uncertainty  in  source  geometry  (Figure  4.20b)  is  usually  larger 
than  the  error  from  the  effects  of  time-dependent  fields  (Figure  4.20a). 

Interestingly,  we  find  source  geometry  errors  at  areas  2  and  8  are  substantially 
larger  than  for  the  other  anomalies.  The  large  error  in  area  2  arises  from  the  existence 
of  two  nearly  equal  local  minima  in  its  DD-CM  error  map,  which  is  not  seen  at  any 
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other  anomaly  (see  Figure  4.18).  The  two  minima  are  approximately  110°  apart. 
Small  displacements  of  the  dipoles  from  their  nominal  locations  flip  the  best-fit 
solution  into  the  other  minimum,  producing  its  large  uncertainty  ellipses,  particularly 
for  day  09096  (Figure  4.20b).  One  of  our  four  best-fit  solutions  (day  99033)  falls  in 
this  secondary  minimum.  Because  area  1  and  2  have  slight  overlap  (Figure  4.1),  we 
tested  if  the  double  minima  could  be  a  result  of  area  1  data  influencing  the  regression 
at  area  2.  We  ran  the  DD-CM  algorithm  against  data  for  area  2,  but  this  time  excluded 
spacecraft  measurements  over  area  1.  The  results  are  not  significantly  different;  the 
two  minima  are  still  seen  for  area  2.  Using  the  GD-VM  algorithm,  one  finds  the  best- 
fit  directions  for  area  2  are  approximately  between  the  two  minima  found  by  the  DD- 
CM  algorithm  (Figure  4.20c),  -60°  from  either  one.  This  difference  between  the  two 
algorithms  is  the  largest  for  any  area.  Further  work  will  be  required  to  better 
understand  the  nature  of  the  error  space  at  area  2,  and  its  source  body  characteristics; 
we  will  investigate,  for  example,  oppositely  magnetized  blocks  of  magnetization  in 
this  region  (cf.  Parker,  1988). 

For  area  8,  the  large  error  arises  largely  from  the  flatness  of  the  error  space, 
instead  of  the  existence  of  multiple  minima  (Figure  4.18).  In  particular,  the 
declination  of  this  nearly  equator-pointing  magnetization  vector  is  poorly  constrained, 
with  results  spanning  a  range  of  -100°  of  arc  (a  result  found  to  be  true  using  either 
the  DD-CM  or  GD-VM  algorithm).  Again,  without  the  error  maps  provided  by  the 
DD-CM  algorithm,  we  would  not  have  been  able  to  isolate  and  understand  these 
sources  of  uncertainty  in  areas  2  and  8. 
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Next,  we  calculate  the  mean  magnetization  strength  for  all  study  areas  using 
the  best-fit  depths  and  magnetic  moments  recovered  by  both  the  DD-CM  and  GD- 
VM  algorithm  (Table  4.1).  For  this  calculation,  the  area  of  the  magnetic  anomalies  is 
taken  as  the  area  across  which  the  magnetic  field  B  >  Bmax/4,  where  Bmax  is  the  peak 
magnetic  field  in  a  given  study  area.  We  choose  to  approximate  the  thickness  of  the 
magnetic  material  as  twice  the  average  of  the  best-fit  magnetic  source  depths.  Most 
values  cluster  in  the  0.1-0.3  A/m  range  found  by  [Purucker  et  al,  2012],  who  studied 
a  region  similar  to  area  1.  A  few  isolated  datasets  at  the  lowest  altitudes  show  higher 
magnetizations,  but  with  the  exception  of  area  9,  average  magnetizations  for  a 
particular  study  area  across  all  datasets  are  close  to  the  0.1 -0.4  A/m  range  suggested 
by  [Wieczorek  and  Weiss,  2010].  Even  for  area  9,  only  the  GD-VM  average 
magnetization  is  outside  this  range. 

In  most  cases  the  GD-VM  algorithm  returns  higher  values  for  magnetizations 
than  the  DD-CM  algorithm.  This  is  a  result  of  the  shallower  depths  returned  by  the 
GD-VM  algorithm,  which  arise  because  of  the  sheet-like  nature  of  the  magnetic 
source,  compared  to  the  dipolar  model  in  the  DD-CM  algorithm.  In  general  the  depth 
of  the  magnetic  source  bodies  is  not  well  constrained  from  our  results,  due  to  the 
possibility  of  non-uniquely  trading  depth  with  moment,  but  ultimately  does  not  factor 
in  to  our  paleo-directional  analysis. 

Finally,  best-fit  paleopoles  for  all  datasets  are  calculated  from  the  best-fit 
magnetization  directions  [Butler,  1998]  (Figure  4.21,  Table  4.3).  The  uncertainties 
derived  from  the  methods  described  in  the  Methods  section  yield  a  circular  error 
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ellipse  for  each  area  that  becomes  distorted  when  its  paleopole  is  calculated.  Hence 
we  calculate  the  two  different  semi-axes  of  the  ellipse  along  the  great  circle  path  from 
the  site  to  pole  {dp)  and  the  semi-axes  of  the  ellipse  perpendicular  to  that  path  {dm) 
(Butler,  1998).  For  area  2,  we  also  show  the  paleopoles  from  the  two  local  error 
minima  found  in  the  DD-CM  method,  particularly  because  one  falls  very  close  to  the 
paleopole  for  area  1  (within  its  uncertainty),  while  the  other  falls  close  to  the 
paleopoles  from  areas  3,  6,  and  8  (within  the  uncertainties  of  6  and  8). 
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Table  4.1.  Lunar  Prospector  (LP)  and  Kaguya  (KG)  datasets  used  for  inversions,  with  best-fit 
magnetic  characteristics  found  using  both  the  DD-CM  and  GD-VM  algorithms.  Mean 
observational  altitude,  data  day  and  site  latitude  and  longitude  are  shown.  GD-VM  dipoles 
column  indicates  the  number  of  gridded  dipoles  used  for  the  GD-VM  simulation.  Calculated 
magnetizations  for  both  algorithms  are  shown.  Highlighted  rows  correspond  to  figures  shown  in 
the  main  text  (Figure  4.8  -  Figure  4.17);  corresponding  GD-VM  dipole  maps  for  these  datasets 
are  also  shown  in  Figure  4.58  -  Figure  4.67.  Best-fit  figures  for  non-highlighted  rows  may  be 
found  in  Figure  4.23  -  Figure  4.57. 


DD-CM  algorithm 

GD-VM  Algorithm 

study  Area 

Spacecraft 

Year 

Day 

Latitude  Longitude 

Alt 

Depth 

Moment 

Inc 

Dec 

Mag 

Dipoles 

Depth 

Inc 

Dec 

Mag 

(deg  N) 

(deg  E) 

(km) 

(km) 

(10^^  Am^)  (deg)  (deg)  (A/m) 

(km) 

(deg) 

(deg)  (A/m) 

1 

LP-MAG 

1999 

142 

-25 

-184 

22.4 

25 

8.4 

13 

-123 

0.12 

580 

5 

18 

-128 

0.70 

LP-MAG 

1999 

61 

-25 

-184 

32.3 

24 

6.4 

32 

-116 

0.09 

580 

5 

22 

-128 

0.62 

K-MAG 

2009 

96 

-25 

-184 

39.4 

36 

11.0 

42 

-132 

0.06 

580 

26 

30 

-126 

0.12 

K-MAG 

2009 

123 

-25 

-184 

33 

37 

6.4 

20 

-114 

0.11 

580 

14 

20 

-116 

0.24 

Merged 

580 

6 

19 

-123 

0.42 

2 

LP-MAG 

1999 

33 

-20 

-188 

32.8 

46 

21.0 

36 

21 

0.06 

1408 

21 

11 

-71 

0.20 

K-MAG 

2009 

69 

-20 

-188 

44.7 

37.5 

26.0 

29 

-129 

0.06 

1408 

39 

5 

-54 

0.10 

K-MAG 

2009 

96 

-20 

-188 

39.4 

40 

15.0 

26 

-135 

0.07 

1408 

38 

7 

-52 

0.11 

K-MAG 

2009 

123 

-20 

-188 

33 

40 

20.0 

32 

-141 

0.07 

1408 

26 

7 

-53 

0.16 

Merged 

1408 

12 

8 

-59 

0.14 

3 

LP-MAG 

1999 

61 

-27 

-192 

32.3 

47 

19.0 

31 

1 

0.08 

432 

21 

20 

3 

0.16 

K-MAG 

2009 

69 

-27 

-192 

44.7 

55 

62.0 

32 

-8 

0.10 

432 

27 

13 

-1 

0.16 

K-MAG 

2009 

96 

-27 

-192 

39.4 

49.5 

18.0 

49 

4 

0.07 

432 

23 

27 

-8 

0.13 

K-MAG 

2009 

123 

-27 

-192 

33 

48 

17.0 

49 

2 

0.06 

432 

23 

21 

2 

0.13 

Merged 

432 

11 

25 

-8 

0.15 

4 

LP-MAG 

1999 

115 

-33 

-183 

32.3 

28 

6.4 

-64 

-46 

0.06 

320 

24 

-55 

-56 

0.09 

LP-MAG 

1999 

142 

-33 

-183 

22.4 

40.5 

12.0 

-73 

-89 

0.07 

320 

13 

-52 

-50 

0.20 

K-MAG 

2009 

96 

-33 

-183 

39.4 

28 

6.2 

-69 

-23 

0.06 

320 

24 

-70 

-28 

0.08 

K-MAG 

2009 

151 

-33 

-183 

13.7 

23.5 

5.7 

-82 

86 

0.06 

320 

8 

-62 

-42 

0.30 

Merged 

320 

9 

-85 

-73 

0.17 

5 

K-MAG 

2009 

96 

-19 

-165 

39.4 

32.5 

11.0 

-33 

-16 

0.06 

850 

24 

-23 

-23 

0.15 

K-MAG 

2009 

151 

-19 

-165 

13.7 

55 

34.0 

-24 

-12 

0.12 

850 

25 

-34 

-33 

0.17 

LP-MAG 

1999 

142 

-19 

-165 

22.4 

46 

21.0 

-29 

-23 

0.09 

850 

26 

-30 

-33 

0.14 

K-MAG 

2009 

123 

-19 

-165 

33 

38.5 

7.9 

-56 

6 

0.06 

850 

28 

-36 

-25 

0.11 

Merged 

850 

12 

-30 

-39 

0.14 

6 

LP-MAG 

1999 

142 

-24 

-165 

22.4 

46 

23.0 

47 

-38 

0.10 

576 

23 

50 

-28 

0.16 

LP-MAG 

1999 

33 

-24 

-165 

32.8 

55 

33.0 

49 

24 

0.08 

576 

31 

39 

53 

0.13 

LP-MAG 

1999 

115 

-24 

-165 

32.3 

55 

30.0 

41 

-8 

0.08 

576 

21 

47 

18 

0.16 

K-MAG 

2009 

96 

-24 

-165 

39.4 

46 

22.0 

49 

31 

0.07 

576 

32 

31 

-21 

0.12 

K-MAG 

2009 

123 

-24 

-165 

33 

55 

21.0 

43 

-14 

0.08 

576 

34 

38 

-25 

0.11 

Merged 

576 

12 

43 

15 

0.14 

7 

LP-MAG 

1999 

33 

-27 

-165 

32.8 

55 

22.0 

35 

-160 

0.09 

341 

26 

28 

-132 

0.15 

LP-MAG 

1999 

115 

-27 

-165 

32.3 

41.5 

12.0 

45 

-180 

0.07 

341 

25 

38 

159 

0.12 

K-MAG 

2009 

96 

-27 

-165 

39.4 

24.5 

8.0 

46 

152 

0.08 

341 

21 

34 

155 

0.14 

K-MAG 

2009 

123 

-27 

-165 

33 

37 

11.0 

40 

131 

0.07 

341 

24 

34 

143 

0.12 

Merged 

341 

14 

39 

165 

0.13 

8 

K-MAG 

2009 

123 

-31 

-164 

33 

25.5 

6.9 

2 

146 

0.05 

1100 

39 

-3 

-131 

0.08 

LP-MAG 

1999 

61 

-31 

-164 

32.3 

46 

21.0 

1 

-172 

0.03 

1100 

45 

3 

134 

0.06 

K-MAG 

2009 

151 

-31 

-164 

13.7 

25 

5.9 

-17 

138 

0.05 

1100 

13 

0 

123 

0.31 

LP-MAG 

1999 

142 

-31 

-164 

22.4 

55 

20.0 

-17 

-128 

0.04 

1100 

36 

-6 

-119 

0.10 

K-MAG 

2009 

96 

-31 

-164 

39.4 

28 

7.8 

-1 

-162 

0.03 

1100 

46 

-3 

-132 

0.07 

Merged 

1100 

13 

-1 

-119 

0.12 

9 

LP-MAG 

1999 

61 

-36 

-198 

32.3 

29.5 

26.0 

1 

119 

0.16 

672 

7 

-4 

145 

0.52 

LP-MAG 

1999 

88 

-36 

-198 

32.2 

34 

37.0 

-12 

146 

0.13 

672 

9 

-11 

162 

0.43 

LP-MAG 

1999 

115 

-36 

-198 

32.3 

25 

30.0 

-6 

139 

0.14 

672 

7 

-8 

141 

0.49 

LP-MAG 

1999 

142 

-36 

-198 

22.4 

22.5 

17.0 

-2 

144 

0.18 

672 

4 

-5 

153 

0.95 

K-MAG 

2009 

96 

-36 

-198 

39.4 

26 

23.0 

-5 

125 

0.16 

672 

7 

-8 

141 

0.41 

K-MAG 

2009 

123 

-36 

-198 

33 

23.5 

15.0 

-5 

135 

0.19 

672 

5 

-3 

143 

0.65 

Merged 

672 

4 

-7 

146 

0.58 

10 

LP-MAG 

1999 

61 

-14 

-153 

32.3 

36 

25.0 

34 

-180 

0.05 

506 

15 

15 

156 

0.20 

LP-MAG 

1999 

88 

-14 

-153 

32.2 

20.5 

8.3 

32 

166 

0.10 

506 

12 

11 

-167 

0.23 

LP-MAG 

1999 

142 

-14 

-153 

22.4 

21 

17.0 

29 

-175 

0.06 

506 

5 

14 

-174 

0.56 

K-MAG 

2009 

123 

-14 

-153 

33 

19 

15.0 

21 

-172 

0.06 

506 

5 

10 

-172 

0.50 

K-MAG 

2009 

151 

-14 

-153 

13.7 

16 

14.0 

25 

-168 

0.07 

506 

4 

10 

-155 

0.87 

Merged 

506 

4.5 

13 

-160 

0.47 
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Table  4.2.  Uncertainty  estimates  for  regression  results.  Minimum  RMS  error  recovered  using  the 
DD-CM  algorithm  for  all  areas  and  altitudes,  the  Fisher  distribution  precision  parameter  A:, 
angular  standard  deviation  (=  81/A:^^^),  obtained  from  Monte  Carlo  simulations  for  the  effects 
of  time-dependent  fields  and  0.5  °  radius  displaced  dipoles.  Altitude-averaged  precision 
parameter  k  and  the  95%  confidence  interval  (0C95)  shown  for  each  area.  Highlighted  rows 
correspond  to  figures  shown  in  the  main  text  (Figure  4.8  -  Figure  4.17). 


Spacecraft 

Year 

Day 

Min.  RMS  error 
(DD-CM) 

(nT) 

Min.  RMS  error 
(GD-VM) 

(nT) 

Noise  simulation 
063  Precision 
parameter  k 

Displaced  dipole  simulation 
063  Precision 

parameter  k 

1 

LP-MAG 

1999 

142 

7.4 

3.9 

1.6 

2526 

lA 

120 

LP-MAG 

1999 

61 

5.6 

2.2 

2.1 

1474 

7.3 

124 

K-MAG 

2009 

96 

2.4 

1.1 

5.0 

264 

4.2 

381 

K-MAG 

2009 

123 

6.3 

2.3 

2.9 

774 

3.8 

446 

2 

LP-MAG 

1999 

33 

5.9 

3.1 

1.5 

3120 

4.5 

332 

K-MAG 

2009 

69 

3.7 

2.3 

3.3 

620 

5.7 

203 

K-MAG 

2009 

96 

4.0 

2.2 

2.1 

1542 

44.1 

3 

K-MAG 

2009 

123 

5.7 

3.2 

1.7 

2318 

14.3 

32 

3 

LP-MAG 

1999 

61 

4.1 

2.0 

3.1 

676 

5.7 

202 

K-MAG 

2009 

69 

3.1 

1.6 

3.1 

674 

7.2 

128 

K-MAG 

2009 

96 

1.7 

0.9 

4.3 

353 

6.1 

177 

K-MAG 

2009 

123 

2.8 

1.3 

5.3 

234 

6.4 

160 

B 

LP-MAG 

1999 

115 

1.5 

1.1 

6.4 

159 

4.7 

300 

LP-MAG 

1999 

142 

4.4 

2.0 

4.2 

375 

5.3 

232 

H 

K-MAG 

2009 

96 

1.1 

0.8 

7.2 

127 

5.5 

219 

K-MAG 

2009 

151 

8.1 

3.7 

1.9 

1761 

8.6 

89 

5 

K-MAG 

2009 

96 

2.0 

1.0 

5.3 

232 

3.7 

472 

K-MAG 

2009 

151 

9.1 

4.0 

1.3 

3959 

3.6 

501 

LP-MAG 

1999 

142 

4.6 

1.9 

2.0 

1648 

3.3 

601 

K-MAG 

2009 

123 

2.2 

1.1 

5.2 

243 

4.1 

383 

6 

LP-MAG 

1999 

142 

5.7 

3.1 

2.4 

1150 

15.1 

29 

LP-MAG 

1999 

33 

3.8 

1.6 

4.8 

288 

11.5 

50 

LP-MAG 

1999 

115 

3.3 

2.0 

4.8 

283 

5.3 

234 

K-MAG 

2009 

96 

2.2 

1.1 

7.8 

109 

6.9 

140 

K-MAG 

2009 

123 

2.7 

1.5 

4.6 

308 

4.4 

343 

7 

LP-MAG 

1999 

33 

3.1 

1.6 

5.0 

262 

4.6 

311 

LP-MAG 

1999 

115 

2.2 

1.3 

6.0 

184 

6.3 

164 

K-MAG 

2009 

96 

1.5 

0.9 

8.9 

83 

8.3 

95 

K-MAG 

2009 

123 

2.4 

1.4 

5.2 

244 

6.7 

147 

8 

K-MAG 

2009 

123 

2.6 

1.7 

4.4 

347 

28.0 

8 

LP-MAG 

1999 

61 

2.5 

1.5 

6.0 

184 

3.9 

443 

K-MAG 

2009 

151 

7.5 

3.9 

2.4 

1121 

36.6 

5 

LP-MAG 

1999 

142 

4.4 

2.6 

3.7 

476 

3.2 

625 

K-MAG 

2009 

96 

1.6 

1.1 

7.2 

128 

5.3 

235 

9 

LP-MAG 

1999 

61 

6.3 

3.2 

1.4 

3496 

4.3 

352 

LP-MAG 

1999 

88 

5.6 

2.3 

1.8 

1999 

8.1 

99 

LP-MAG 

1999 

115 

5.7 

2.1 

1.9 

1837 

9.6 

71 

LP-MAG 

1999 

142 

9.0 

3.5 

1.0 

6997 

8.5 

91 

K-MAG 

2009 

96 

4.0 

1.5 

2.8 

853 

8.1 

101 

K-MAG 

2009 

123 

4.5 

1.5 

1.8 

2004 

5.4 

228 

10 

LP-MAG 

1999 

61 

3.0 

1.5 

3.0 

730 

7.7 

112 

LP-MAG 

1999 

88 

3.2 

1.7 

7.9 

104 

11.9 

46 

LP-MAG 

1999 

142 

4.7 

2.3 

2.2 

1315 

10.2 

63 

K-MAG 

2009 

123 

1.9 

1.3 

3.9 

431 

7.3 

124 

K-MAG 

2009 

151 

8.0 

4.2 

1.2 

4571 

15.5 

27 

86 


Table  4.3.  Errors  and  paleopoles  arising  from  using  variable  altitude  data.,  The  paleolatitude, 
paleolongitude,  and  value  of  k  and  the  95%  confidence  ellipse  (0C95)  obtained  by  combining 
the  best-fit  directions  of  the  DD-CM  and  GD-VM  algorithm  best-fit  directions.  The  0C95  value  is 
added  to  the  mean  dispersions  in  Figure  20a  and  20b  to  obtain  the  final  circular  error,  also 
shown.  ,The  distorted  version  of  this  error  circle  at  the  paleopole  is  also  shown  (Figures  21  and 
22):  the  semi-axis  of  the  ellipse  along  the  great  circle  path  from  site  to  pole  (dp)  and  semi-axis  of 
the  ellipse  perpendicular  to  that  path  (dm). 


Area 

Altitude  averaged 

k  a-95 

Paleolatitude 

Paleolongitude 

Final 

circular  error 

dp 

dm 

58.5 

6.8 

-34.8 

84.2 

12.2 

7.0 

13.1 

3.3 

33.9 

8.5 

-259.6 

41.0 

24.0 

44.3 

38.7 

8.4 

46.9 

-194.9 

15.7 

9.6 

17.4 

27.8 

9.9 

-51.8 

-140.9 

16.4 

24.5 

28.4 

33.0 

9.1 

-68.1 

-75.5 

13.3 

8.6 

15.1 

14.5 

12.4 

38.0 

-165.1 

19.9 

16.6 

25.7 

7 

11.2 

16.1 

-78.5 

-62.0 

24.5 

18.2 

29.9 

8 

3.2 

30.9 

54.2 

-142.2 

37.1 

18.7 

37.3 

B 

50.0 

5.9 

37.0 

-247.6 

12.6 

6.3 

12.6 

25.4 

9.2 

-83.3 

82.9 

18.9 

10.3 

19.8 

4.4.2  Comparison  with  other  magnetic  paleopoles 


In  Figure  4.22b-d  we  show  paleopoles  from  isolated  magnetic  anomalies 
studied  by  Takahashi  et  al.  (2014)  and  Arkani-Hamed  and  Boutin  (2014),  as  well  as 
paleopoles  inferred  from  the  magnetization  of  the  crust  at  craters  [Arkani-Hamed  and 
Boutin,  2014],  The  magnetized  crust  paleopoles  from  Arkani-Hamed  and  Boutin 
(2014)  are  the  means  of  the  values  listed  in  their  Table  2.  None  of  the  anomalies 
studied  by  these  groups  were  inside  SPA.  To  compare  our  results  with  these  three 
datasets,  we  reverse  all  of  our  paleopoles  into  the  same  southern  hemisphere  (Figure 
4.22a),  as  in  Takahashi  et  al.  (2014).  Here  we  have  dropped  the  paleopole  found 
from  the  mean  of  all  data  sets  at  area  2,  and  only  show  the  paleopoles  from  the  two 
local  error  minima  in  the  DD-CM  algorithm. 
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Takahashi  et  al.  (2014)  found  paleopoles  that  cluster  into  two  groups,  one  near 
the  present  pole,  and  another  at  mid-latitudes.  We  find  a  wider  dispersion  in  the 
distribution  of  our  paleopoles,  but  we  do  find  some  clustering  at  the  present  pole 
(areas  5,  7  and  10),  as  well  some  near  their  mid-latitude  cluster  (areas  3,  6,  and  8). 
The  remaining  poles  are  not  easily  assigned  to  either  of  these  clusters.  However, 
overall,  the  paleopoles  we  find  do  seem  to  avoid  longitudes  on  the  farside  (the  bottom 
half  of  the  sphere  in  Figure  4.22a),  and  latitudes  <  30°.  None  of  the  paleopoles  from 
Arkani-Hamed  and  Boutin  (2014)  show  obvious  correlation  with  any  of  the  clustering 
found  here  or  in  Takahashi  et  al.  (2014). 

Finally,  we  note  that  many  of  the  paleopoles  are  substantially  separated  by 
their  error  ellipses,  such  that  it  is  unlikely  that  only  one  or  two  paleopoles  would  be 
consistent  with  all  of  the  anomalies,  without  severely  affecting  their  RMS  error 
values.  This  can  also  be  visualized  by  examining  the  error  spaces  in  Figure  18  and  19; 
because  most  anomalies  are  at  similar  latitudes  and  longitudes  (most  within  -30°),  the 
diversity  in  paleopole  locations  is  mostly  determined  by  the  diversity  in 
magnetization  directions.  Hence,  the  diversity  in  RMS  error  minima  in  Figure  18/19 
also  graphically  illustrates  the  diversity  in  paleopole  locations. 

In  the  next  section,  we  address  possible  origins  for  the  dispersion  seen  in  our 
study,  and  test  some  of  the  hypotheses  for  anomaly  formation. 
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4.5  Discussion 


4.5.1  Paleopole  directions 

The  wide  variation  in  magnetic  paleopoles  derived  from  SPA’s  magnetic 
anomalies  presents  a  puzzle.  In  Sections  4.5.1  -  4.5.4,  we  discuss  the  implications  of 
the  diverse  paleopole  locations  for  the  formation  of  these  anomalies,  the  history  of 
true  polar  wander,  and  the  nature  of  the  lunar  dynamo. 

4.5.2  Magnetization  by  Impact  Processes 

A  long-standing  hypothesis  is  that  the  strongest  lunar  magnetic  anomalies  are 
genetically  related  to  the  antipodes  of  the  Imbrium,  Orientale,  Serenitatis  and  Crisium 
basins  [Moore  et  ah,  1974;  Hood  and  Williams,  1989;  Lin  et  al,  1998;  Hood  et  al, 
2001].  Compression  and  amplification  of  the  interplanetary  magnetic  field  (IMF)  by 
impact-produced  plasma  may  be  strongest  at  the  basin  antipode,  where  impact  ejecta 
may  also  preferentially  collect  [Hood  and  Artemieva,  2008].  The  antipodes  to  the 
Imbrium  and  Serenitatis  basins  are  close  to  our  study  areas  (Figure  4.1),  suggesting 
this  process  may  be  responsible  for  forming  the  anomalies  we  have  examined  here.  If 
true,  antipodal  ejecta  may  become  magnetized  via  either  thermo-remanent 
magnetization  (TRM)  or  shock-remanent  magnetization  (SRM).  Below  we  assess 
these  two  possibilities. 


89 


4.5.2. 1  TRM  in  an  impact-produced  field 

If  the  magnetization  were  produced  by  a  TRM,  hot  ejecta  would  cool  in  the 
presence  of  transient  IMF-amplified  fields  that  would  last,  at  most,  for  one  day. 
Assuming  a  thermal  diffusivity  of  10'  m  /s,  the  thermal  cooling  length-scale  for  one 
day  is  ~1  meter.  Therefore,  TRM  would  be  restricted  to  1  meter  of  material.  The 
magnetic  moments  we  find  across  all  SPA  study  areas  range  between  lO'^  Am^  and 
10  Am  .  Using  source  body  horizontal  extents  from  Figure  4.1  (black  dashed 
boxes),  depths  of  ~  1  m  lead  to  high  TRM  magnetizations  ranging  between  10^  to  lO"^ 
A/m.  These  values  are  3-4  orders  of  magnitude  higher  than  samples  recovered  by 
Apollo  missions.  Additionally,  the  top  ~1  m  would  have  been  completely  overturned 
and  demagnetized  in  the  time  since  the  antipodal  impact  ~4  billion  years  ago  [Arnold, 
1975].  Therefore,  TRM  from  impact-related  fields  is  implausible,  in  agreement  with 
[Hood  and  Artemieva,  2008]. 

4.5.2.2  SRM  in  an  impact-produced  field 

Alternately,  impact  shock  pressures  from  the  deposition  of  ejecta  at  the 
antipode  may  create  SRM  in  the  ejecta  deposit.  However,  to  allow  unidirectional 
SRM  to  be  stably  imparted  to  a  large-scale  (e.g.  >  4000  km^  for  Areas  1  and  2; 
>30000  km^  for  Areas  5-8)  geologic  formation,  the  shock  waves  must  pass  through 
the  rock  when  the  entire  region  is  at  rest.  Shock  waves  begin  propagating  through  the 
rock  instantly  upon  impact,  yet  at  that  moment  the  ejecta  is  still  traveling  at  2.0-2.4 
km/s  [Hood  and  Artemieva,  2008].  Therefore,  the  ejecta  at  the  antipode  will  change 
orientation  as  it  comes  to  rest,  after  any  SRM  has  formed,  such  that  the  total 
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remanence  of  the  ejecta  deposit  will  be  randomized  and  nulled.  Hood  and  Artemieva 
(2008)  also  rule  out  SRM  in  the  ejecta,  using  a  different  line  of  reasoning. 

It  is  possible  that  the  underlying  rock  could  be  shocked  by  the  impacting 
ejecta  and  acquire  a  unidirectional  SRM  [Hood  and  Artemieva,  2008].  However,  we 
suggest  that  there  are  two  problems  with  this  hypothesis.  Firstly,  one  would  expect 
nearly  vertical  local  magnetization  if  the  SRM  producing  field  was  compressed  by 
solar  wind  plasma  (see,  for  example,  field  lines  in  Figure  8-9,  Hood  and  Artemieva, 
2008).  Instead,  only  one  magnetic  anomaly,  area  4,  has  a  nearly  vertical 
magnetization,  and  the  rest  show  a  preference  for  low  inclinations,  if  anything  (the 
rest  are  at  least  >33°  away  from  the  vertical.  Figure  4.20). 

Secondly,  deposition  of  impact  ejecta  would  take  place  over  a  short  period  of 
time,  over  which  the  ambient  field  direction  is  likely  to  be  nearly  constant.  There  are 
only  two  basin  antipodes  within  SPA,  which  implies  that  the  magnetic  anomalies 
should  have  one  of  two  magnetization  directions.  Roughly,  areas  1-3  (northwest 
cluster),  4,  and  9  are  antipodal  to  Imbrium,  while  areas  5-8  (eastern  stripes),  and  10 
are  antipodal  to  Serenitatis.  First,  we  consider  the  five  anomalies  at  the  Imbrium 
antipode.  Here  we  find  that  areas  I  and  3  have  magnetization  vectors  109°  apart.  The 
magnetization  at  area  4  is  85°  from  that  at  area  I  and  109°  from  area  3.  The 
magnetization  at  area  9  (Mare  Ingenii)  is  76°  from  area  I,  40°  from  area  3,  and  128° 
from  area  4.  All  of  these  separations  are  well  outside  the  error  ellipses.  The  only  two 
clusters  are  the  directions  of  areas  1  and  2,  and  possibly  area  3  with  the  second  error 
minimum  obtained  at  area  2  (here  we  do  not  include  the  extent  of  the  error  ellipse  of 
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area  2,  which  just  overlaps  with  area  9,  since  that  error  ellipse  is  affected  by  the 
existence  of  two  error  minima).  In  sum,  there  are  at  least  four  widely  separated 
magnetization  directions  near  the  Imbrium  antipode,  all  well  separated  by  their  error 
ellipses.  It  is  plausible  the  similar  directions  at  areas  1  and  2  suggest  this  pair  was 
magnetized  by  deposition  of  impact  ejecta,  but  it  would  not  explain  the  magnetization 
of  the  other  areas. 

At  the  five  anomalies  at  the  Serenitatis  antipode,  areas  5,  7,  and  10  have 
similar  directions,  as  do  areas  6  and  8,  but  the  two  groups  are  over  -120°  apart,  and 
well  separated  by  error  ellipses.  Therefore,  we  conclude  again  that  it  is  not  likely  that 
the  Serenitatis  impact  is  responsible  for  magnetizing  all  of  these  anomalies. 

In  sum,  the  diversity  of  directions  argues  against  the  SRM  hypothesis,  or  at 
least  allows  only  a  subset  of  the  geographically  clustered  anomalies  to  be  due  to 
SRM.  Similar  arguments  can  be  applied  to  ruling  out  SRM  from  surface  seismic 
waves.  Considering  all  of  the  observations  above,  we  conclude  that  at  least  some  of 
the  magnetic  source  bodies  in  SPA  were  magnetized  in  a  lunar  dynamo  field,  rather 
than  a  field  associated  with  impact  events. 

4.5.3  Magnetized  South  Pole-Aitken  basin  ejecta  or  volcanic 
bodies? 

[Wieczorek  et  al,  2012]  proposed  that  material  from  the  SPA  impactor  impact 
might  be  the  source  of  many  of  the  SPA  basin’s  magnetic  anomalies,  and  even  other 
anomalies  across  the  Moon.  Under  this  hypothesis,  hot  iron-enriched  material  from 
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the  SPA  impactor  acquired  a  TRM  in  a  dynamo  field.  However,  if  SPA  ejecta  landed 
hot  and  cooled  in  a  dynamo  field,  the  resulting  anomalies  should  all  have  the  same 
magnetization  direction  (but  see  section  4.5. 3.4).  Instead,  the  diverse  magnetization 
directions  (see  section  4.5. 1.2,  above)  suggest  that  they  were  magnetized  at  different 
epochs. 

The  last  remaining  viable  hypothesis  is  that  magnetic  anomalies  in  the 
northern  SPA  basin  formed  as  a  result  of  magnetized  sub-surface  dikes  [Purucker  et 
al,  2012].  Dikes  forming  over  long  time  periods  would  permit  different 
magnetization  directions  during  different  magnetic  epochs.  The  cluster  of  three 
paleopoles  close  to  the  present  pole  (areas  5,  7  and  10,  Figure  4.22a)  would  suggest  a 
traditional  axial-aligned  dynamo  magnetized  these  dikes  when  the  Moon  was  in  its 
present  orientation.  However,  the  diversity  in  paleopoles  seen  is  still  enigmatic  (see 
Section  4.5.3).  Further,  if  the  dike  hypothesis  is  true,  this  implies  that  the  dikes  near 
SPA  are  special  in  some  way,  since  the  nearside  of  the  Moon,  covered  much  more 
extensively  by  volcanism  (and  presumably  associated  with  subsurface  dikes),  shows 
no  magnetic  structures  like  those  at  SPA.  Andrews-Hanna  et  al.  (2014)  reported  linear 
gravity-gradient  anomalies  that  may  be  dikes,  but  so  far,  no  obvious  correlation 
between  these  structures  and  magnetic  anomalies  have  been  found.  We  note  that 
Gong  and  Wieczorek  (2016)  find  a  correlation  between  magnetization  and  gravity 
anomalies  in  some  locations. 
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4.5.4  Possible  Explanations  for  Diverse  Magnetic  Field 
Directions 

The  above  arguments  suggest  the  magnetizing  fields  for  many  of  the  SPA 
anomalies  arose  from  a  dynamo.  However,  they  do  not  offer  an  obvious  explanation 
for  the  diversity  in  field  directions.  Due  to  the  small  size  of  the  lunar  core,  and  the 
rapid  decay  of  magnetic  quadrupole  and  higher  terms  as  a  function  of  distance,  it  is 
likely  that  the  ancient  lunar  dynamo  was  dominantly  dipolar  at  the  surface  [Weiss  and 
Tikoo,  2014].  If  the  dipole  was  aligned  with  the  Moon’s  spin  axis,  the  magnetization 
directions  contain  information  about  the  Moon’s  paleopole.  The  diverse  paleopoles  in 
Figure  4.21  seem  to  imply  large  amounts  of  true  polar  wander  [Goldreich  and 
Toomre,  1969;  Runcorn,  1983].  While  the  diversity  in  Figure  4.21  is  surprising,  in 
comparison  with  the  results  of  [Arkani-Hamed  and  Boutin,  2014;  Takahashi  et  al, 
2014]  (Figure  4.22a),  the  diversity  of  paleopoles  seen  in  our  findings  is  still  greater. 
Below  we  consider  some  possible  explanations  for  the  diverse  paleopole  locations. 

4.5.4.1  Non-axially  aligned  dipoles  and  impact-induced  dynamos 

Currently,  our  understanding  of  the  nature  of  the  lunar  dynamo  is  limited. 
Paleomagnetic  studies  favor  a  dynamo  that  existed  from  approximately  4.2-3. 6  Ga 
ago.  Mechanisms  for  sustaining  such  a  long-lived  core  dynamo  are  uncertain,  with 
recent  proposals  for  dynamos  driven  by  mechanical  stirring  from  impacts  and 
precession  [Dwyer  et  ah,  2011;  Le  Bars  et  ah,  2011;  Weiss  and  Tikoo,  2014]. 
However,  it  is  not  known  if  these  dynamos  produce  the  same  field  organization  as  the 
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Earth’s  dynamo.  It  may  be  possible  that  these  exotic  dynamos  exhibit  more  variable 
dipole  axis  directions,  which  could  explain  the  diversity  of  paleopoles  seen. 
Unfortunately,  more  work  modeling  small  dynamos  like  that  of  the  Moon  is  needed 
before  evaluating  this  hypothesis  further. 

4.5.4.2  Secular  Variation 

Another  possible  origin  for  the  diverse  paleopole  locations  is  geomagnetic 
secular  variation  of  a  dynamo  that  is  on  average  aligned  with  the  lunar  spin  axis.  The 
Earth’s  spin  axis  is  presently  11°  away  from  its  magnetic  dipole  axis,  and  it  is 
plausible  that  the  Moon’s  ancient  dipole  axis  was  also  not  exactly  aligned  with  its 
spin  axis.  It  is  also  possible  that  secular  drift  might  be  larger  on  a  body  with  a  small 
core  such  as  the  Moon.  However,  some  of  the  paleopole  locations  at  SPA  are  over 
-75°  apart  (after  accounting  for  the  possibility  of  reversals,  Eigure  4.22a),  and  over 
the  timescale  of  magnetic  anomaly  formation,  the  mean  dipole  orientation  might 
average  to  be  aligned  with  the  spin  axis.  Presently,  we  have  little  information 
available  about  the  lunar  dynamo  to  further  evaluate  the  role  of  secular  variation  in 
explaining  the  diversity  of  paleopole  directions. 

4.5.4.3  True  Polar  Wander  over  long  timescales 

The  large  amount  of  polar  wander  implied  by  the  paleopoles  is  difficult  to 
reconcile  with  other  geophysical  constraints  for  the  orientation  history  of  the  Moon. 
Currently  there  are  two  comprehensive  studies  that  derive  estimates  of  the  degree  of 
polar  wander  on  the  Moon  due  to  long-term  changes  in  the  Moon’s  moments  of 
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inertia  \Garrick-Bethell,  2016].  The  first  uses  the  shape  and  gravity  of  the  Moon 
exterior  to  large  basins  to  establish  the  earliest  orientation  of  the  Moon  [Garrick- 
Bethell  et  al,  2014b].  The  authors  find  the  present  lunar  pole  has  changed  by  -36° 
from  its  earliest  axis,  but  it  does  not  coincide  with  any  of  the  paleopole  clusters  we 
find  (Figure  4.22a).  The  second  study  uses  polar  hydrogen  deposits  as  a  constraint  on 
the  history  of  polar  wander,  and  infers  that  up  to  -10°  of  paleolatitude  change  has 
occurred  [Siegler  et  al,  2016].  They  infer  large  longitude  changes,  but  they  cannot 
produce  the  large  paleopole  changes  implied  by  the  anomalies  studied  here.  In 
summary,  the  limited  number  of  available  geophysical  models  that  estimate  lunar 
polar  wander  cannot  produce  the  diversity  and  large  magnitude  of  paleopole  changes 
required  to  explain  our  observations.  Interestingly,  we  find  that  all  of  our  paleopoles 
are  >35°  from  the  present  equator  (Figure  22a,  accounting  for  reversals),  which  may 
be  due  to  the  difficulty  in  producing  the  large  required  changes  in  the  Moon’s 
moments  of  inertia. 

4.5.4.4  True  Polar  Wander  due  to  SPA  formation 

Very  large  initial  changes  in  the  Moon’s  moments  of  inertia  due  to  SPA’s 
crater  might  have  produced  large  amounts  of  polar  wander.  Eventually,  these  changes 
must  have  subsided  over  millions  to  billions  of  years,  since  the  gravity  signature 
observed  today  is  muted  \Zuber  et  al,  2013].  If  hot  material  was  cooling  in  the 
presence  of  a  dynamo  field  throughout  the  Moon’s  reorientation,  its  magnetization 
could  in  principle  capture  multiple  lunar  orientations,  thereby  producing  the  diversity 
of  paleopole  locations.  Some  of  this  hot  material  could  be  iron-rich  material  from  the 
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SPA  impactor  [Wieczorek  et  al,  2012],  Following  impact,  the  free  precession 
damping  time  of  the  Moon  is  -2x10^  years  \Peale,  1976;  Williams  et  al,  2001], 
assuming  a  semimajor  axis  of  30  Earth  radii,  dissipation  quality  factor  Q  =  5Q  and 
degree-2  love  number  ^2  =  0.1  (which  we  assume  to  be  representative  of  the  Moon 
when  SPA  formed).  The  length  scale  for  cooling  in  -2x10^  years  is  ~3  km  (assuming 
a  thermal  diffusivity  of  10'^),  which  is  small  compared  to  the  scale  of  the  anomalies 
we  observe,  but  perhaps  not  so  much  as  to  preclude  recording  a  measurable 
magnetization. 

However,  using  a  simple  model  for  reorientation,  we  find  that  density 
anomalies  produced  by  SPA’s  crater  do  not  produce  paleopoles  that  overlap  with 
those  of  its  magnetic  anomalies.  In  our  model,  we  replace  SPA’s  gravity  potential 
inside  the  outer  topography  rim  with  values  between  -3  times  the  maximum,  and  +2 
times  the  maximum  of  the  present  day  potential,  in  increments  of  0.25  times  the 
maximum  potential.  The  inertia  tensor  and  paleopoles  were  then  calculated  for  each 
of  these  cases  using  the  resulting  globally  calculated  degree-2  spherical  harmonics 
(cf.  Garrick-Bethell  et  al.,  2014).  These  positive  and  negative  gravity  values  represent 
very  different  models  of  SPA’s  effects,  but  illustrate  the  range  of  paleopoles  that  are 
possible.  We  find  the  maximum  extent  of  true  polar  wander  from  SPA’s  formation 
passes  close  to  the  paleopoles  from  areas  7  and  10  (Figure  4.22a,  the  negative 
anomaly  path  approaches  the  center  of  SPA  as  the  anomaly  size  grows,  as  expected). 
However,  these  anomalies  are  already  close  to  the  present  pole,  and  do  not  need  to  be 
explained  by  polar  wander.  We  do  find  that  the  extreme  limits  of  polar  wander 
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approach  the  error  ellipses  of  areas  5  and  8,  such  that  this  process  could  plausibly 
explain  some  of  these  magnetic  paleopoles.  Unfortunately,  the  paleopoles  from  areas 
1,  2,  3,  5,  and  9,  and  to  a  lesser  extent  area  6,  do  not  come  close  to  the  paleopoles 
produced  by  SPA.  Thus,  cooling  of  SPA  ejecta  deposits  over  short  timescales  (~10^ 
years),  or  even  dikes  that  formed  when  the  Moon  resided  at  any  of  the  SPA-produced 
paleopoles,  cannot  fully  explain  the  observed  diversity  of  magnetic  paleopoles  we 
find. 

There  are  many  unknowns  in  modeling  the  paleopoles  allowed  by  SPA’s 
formation.  For  example,  reorientation  and  magnetic  anomaly  formation  depends  on 
the  SPA  impact’s  effect  on  dynamo  operation  [Arkani-Hamed  and  Olson,  2010a, 
2010b],  the  inertia  tensor  of  the  Moon  just  before  SPA  formation  and  the  length  scale 
(and  thereby  cooling  timescale)  of  the  materials  making  up  the  anomalies. 
Furthermore,  the  spin  vector  of  the  Moon  will  be  freely  precessing  around  its  angular 
momentum  vector  during  over  the  damping  timescale  of  reorientation.  If  the  dynamo 
after  SPA  impact  retained  its  alignment  along  the  Moon’s  angular  momentum  vector, 
this  precession  could  broaden  the  range  of  paleopoles  permitted  (essentially  accessing 
a  range  of  paleopoles  around  the  paleopoles  shown  in  Figure  4.22a,  with  the  range 
depending  on  the  precession  angle).  Therefore,  we  cannot  definitively  rule  out  polar 
wander  processes  as  the  origin  for  some  of  the  diverse  paleopole  locations  at  SPA. 

A  variant  of  the  hypothesis  of  iron-rich  SPA  ejecta  [Wieczorek  et  al,  2012]  is 
the  cooling  of  iron-rich  “sesquinary”  [Zahnle  et  al.,  2008]  impactors  formed  by  the 
SPA  impact.  Ejected  into  orbit  immediately  after  impact,  studies  on  various  planetary 
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bodies  show  that  these  impactors  can  return  to  a  body  from  10“^  years  {Nayak  et  al, 
2016b]  to  10^  years  [Gladman  et  al.,  1995]  post-impact,  at  approximately  escape 
velocity  [Nayak  and  Asphaug,  2016].  With  reimpacts  spread  across  this  timescale,  we 
propose  that  reaccreting  iron-rich  material  originally  from  SPA,  either  still  hot  from 
the  impact  or  heating  upon  re-impact,  could  record  a  diverse  set  of  orientations  as  the 
Moon  reorients  in  response  to  moment  of  inertia  changes.  Across  this  timescale, 
impact  locations  become  random  [Nayak  et  al,  2016a];  such  sesquinary  magnetism 
caused  by  these  iron-rich  impactors  would  be  widely  distributed  around  the  Moon. 

4.6  Conclusion 

Using  two  different  inversion  methods,  we  find  diverse  directions  of 
magnetization  among  magnetic  anomalies  in  the  northern  SPA  basin.  The  diverse 
directions  help  rule  out  impact-related  fields  as  their  only  origin.  Intrusive  bodies, 
such  as  the  dikes  proposed  by  [Purucker  et  al,  2012],  are  a  plausible  explanation. 
The  diverse  paleopole  locations  could  imply  large  amounts  of  true  polar  wander,  but 
true  polar  wander  inferred  independently  from  gravity  [Garrick-Bethell  etal,  2014b] 
and  hydrogen  deposits  [Siegler  et  al,  2016]  implies  more  modest  changes  in  the 
Moon’s  orientation.  The  diverse  directions  argue  against  the  hypothesis  that  they  all 
formed  from  iron-rich  SPA  ejecta  that  cooled  in  a  dynamo  field  [Wieczorek  et  al, 
2012].  A  simple  gravity  anomaly  model  for  large  amounts  of  SPA-produced 
reorientation  fails  to  explain  at  least  five  of  the  paleopoles  observed,  but  many 
unknowns  remain  in  modeling  this  process.  Some  SPA  ejecta  may  have  produced 
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iron-rich  sesquinary  impactors  that  landed  across  the  Moon  and  recorded  orientation 
changes  as  they  cooled  in  a  dynamo  field.  A  dynamo  that  was  not  aligned  with  the 
lunar  spin  axis  remains  a  plausible  explanation  for  all  observations,  but  gaps  remain 
in  our  ability  confirm  this  hypothesis.  The  wide  variety  of  viable  hypotheses  and  large 
number  of  unknowns  highlight  the  complexity  of  interpreting  the  origins  of  lunar 
magnetic  anomalies  and  their  paleopoles. 
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4.7  Figures 
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Figure  4.1.  Lunar  South  Pole-Aitken  basin  study  areas.  The  SPA  rim  is  per  Garrick-Bethell  and 
Zuber  (2009).  (a)  Magnetic  contours  superimposed  over  LOLA  topography  data  (contours  are 
taken  from  the  map  in  part  b).  Contours  shown  are  5,  8  and  11  nT;  (b)  Magnetic  field  map  with 
study  areas  highlighted  (dashed  black  lines).  Magnetic  field  data  are  from  the  Kaguya  spacecraft 
magnetometer  measurements  taken  on  day  96  in  2009;  at  a  mean  altitude  of  39.4  km.  Black 
crosses  indicate  the  location  of  the  Imbrium  (west)  and  Serenitatis  (east)  basin  antipodes,  (c) 
Magnetic  field  map  taken  from  a  spherical  harmonic  model  [Purucker  and  Nicholas,  2010]  for  the 


101 


SPA  region.  Study  areas  are  highlighted.  We  refer  to  areas  1-3  as  the  northwestern  cluster,  and 
areas  5-8  as  the  eastern  stripes. 
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Figure  4.2.  Overview  of  all  datasets  used  for  SPA  areas  1-5.  Each  row  denotes 
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Figure  4.3.  Overview  of  all  datasets  used  for  SPA  areas  6-10.  Each  row  denotes  a  study  area  from 
Figure  4.1.  Mean  measurement  altitudes  are  shown.  Data  collection  days  may  be  read  in  the 
format  YY-DDD,  e.g.,  09123  is  2009  day  123,  and  so  on. 
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Figure  4.4.  Magnetic  field  maps  and  histograms  of  field  strengths  for  the  Mare  Imbrium  region 
in  the  lunar  wake.  Panels  (a)  and  (c)  are  magnetic  field  maps  of  Imbrium  representing  data  from 
LP-MAG  observation  days  118  and  131  in  1999  respectively.  White  dots  indicate  observation 
locations.  The  mean  observation  altitude  is  31.3  km  for  day  118  and  29.8  km  for  day  131.  Panels 
(b)  and  (d)  are  corresponding  histograms  of  the  strength  of  the  total  magnetic  field. 


Figure  4.5.  Error  estimation  due  to  time  variable  (non-crustal)  magnetic  fields.  Kaguya  magnetic 
field  data  for  area  1,  2009  day  96  (left  column),  the  effect  of  random  noise  with  a  standard 
deviation  of  1  nT  added  to  that  data  (middle  column)  and  one  of  our  Monte  Carlo  regressions  for 


104 


the  depth,  moment  and  direction  values  for  data  +  noise  with  the  DD-CM  algorithm  (right 
column).  Dipole  source  geometry  is  identical  to  the  placement  in  Figure  4.8  (black  dots). 


Figure  4.6.  RMS  error  map  for  the  magnetization  direction  for  area  1,  2009  day  96,  centered  at 
0°  inclination  and  0°  declination  (whole-sphere  Mollweide  projection,  in  which  the  southern 
hemisphere  is  positive  inclination).  The  small  white  circular  contour  indicates  one  angular 
standard  deviation  of  dispersion,  from  Monte  Carlo  simulations  of  the  addition  of  1  nT  Gaussian 
noise  to  observations  (Figure  4.5).  The  larger,  outer  white  dashed  contour  indicates  the  minimum 
error  solution  plus  1  nT. 


Figure  4.7.  Error  estimation  due  to  uncertainty  in  source  body  geometry.  Black  points  are  the 
nominal  dipole  locations  for  a  DD-CM  regression  at  area  1.  Circles  represent  the  locations  of 
dipoles  placed  randomly  in  our  Monte  Carlo  error  simulations  and  are  0.5°  in  radius.  K-MAG 
total  field  observations  in  the  lunar  wake  on  2009  day  96.  Compare  to  Figure  4.8. 
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Figure  4.8.  Area  1  best  fit  results  for  both  inversion  algorithms.  (Left  column)  K-MAG 
observations  in  the  lunar  wake  on  2009  day  96,  compared  to  (Middle  column)  a  grid-based  model 
with  variable  magnetizations  (GD-VM)  and  (Right  column)  a  model  with  manually  placed 
dipoles  of  equal  magnetization  (DD-CM).  Source  dipoles  in  the  DD-CM  algorithm  are 
approximately  located  at  maxima  in  the  observed  total  field  (black  dots).  White  dots  in  the  east, 
north  and  radial  panels  are  locations  of  spacecraft  magnetometer  measurements.  The  mean 
measurement  altitude  is  shown  in  Figure  4.2. 
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Figure  4.9.  Area  2  best  fit  results  with  both  inversion  algorithms.  Figure  details  are  as  in  Figure 
4.8,  except  that  magnetometer  measurements  are  taken  from  K-MAG,  2009  day  123. 
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Figure  4.10.  Area  3  best  fit  results  with  both  inversion  algorithms.  Figure  details  as  in  Figure  4.9. 
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Figure  4.11.  Area  4  best  fit  results  with  both  inversion  algorithms.  Figure  details  as  in  Figure  4.8. 
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Figure  4.12.  Area  5  best  fit  results  with  both  inversion  algorithms.  Figure  details  as  in  Figure  4.9. 
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Figure  4.13.  Area  6  best  fit  results  with  both  inversion  algorithms.  Figure  details  as  in  Figure  4.9 
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Figure  4.14.  Area  7  best  fit  results  with  both  inversion  algorithms.  Figure  details  as  in  Figure  4.9 
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Figure  4.15.  Area  8  best  fit  results  with  both  inversion  algorithms.  Figure  details  as  in  Figure  4.9. 
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Figure  4.16.  Area  9  best  fit  results  with  both  inversion  algorithms.  Figure  details  as  in  Figure  4.9. 
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Figure  4.17.  Area  10  best  fit  results  with  both  inversion  algorithms.  Figure  details  as  in  Figure 
4.9. 
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Figure  4.18.  RMS  error  map  for  the  magnetization  direction  for  all  datasets  for  SPA  study  areas 
1-5,  obtained  using  the  DD-CM  algorithm.  Southern  hemisphere  is  positive  inclination;  the 
whole-sphere  Mollweide  projection  is  centered  at  0°  inclination  and  0°  declination.  The  solid 
white  circular  contour  indicates  one  angular  standard  deviation  of  dispersion  from  the  best-fit 
solution  (red  star,  Table  4.1),  from  Monte  Carlo  simulations  of  the  addition  of  1  nT  Gaussian 
noise  to  observations  in  the  DD-CM  algorithm.  The  larger,  dashed  white  contour  indicates  the 
minimum  error  solution  plus  1  nT,  a  measure  of  uncertainty  using  an  arbitrary  threshold 
defined  by  the  measurement  uncertainty  (not  used  in  our  final  analysis).  The  error  from  Monte 
Carlo  simulations  of  the  effects  displacing  the  nominal  dipoles  by  0.5°  is  not  shown,  and  nor  is 
the  95%  confidence  interval  from  using  different  altitude  data  sets. 
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Figure  4.19.  RMS  error  map  for  the  magnetization  direction  for  all  datasets  for  SPA  study  areas 
6-10,  obtained  using  the  DD-CM  algorithm.  Description  is  as  in  Figure  4.18. 
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Figure  4.20.  Best-fit  magnetization  directions  for  study  areas  defined  in  Figure  1,  for  all  datasets 
in  Figure  4.2  and  Figure  4.3.  Directions  derived  from  different  altitudes  and  different  spacecraft 
at  the  same  study  area  are  represented  by  the  same  color.  Positive  inclinations  are  in  the  lower 
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hemisphere  (filled  circles).  (A  &  B)  Results  for  the  DD-CM  algorithm.  Ellipses  in  A  represent 
one  standard  deviation  of  dispersion  from  1  nT  noise  simulations.  Ellipses  in  B  represent  one 
standard  deviation  of  dispersion  from  simulations  of  displaced  dipoles  (0.5°),  (C)  Results  from 
the  GD-VM  algorithm,  including  results  from  regressions  to  merging  data  at  all  altitudes  shown 
in  Figure  4.2  and  Figure  4.3,  within  a  given  study  area  (stars). 
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Figure  4.21.  SPA  magnetic  paleopoles.  Paleopoles  are  the  mean  of  the  combined  results  from  the 
DD-CM  and  GD-VM  algorithms,  at  all  directions  (altitudes)  shown  in  Figure  4.20.  Ellipses 
represent  the  95%  confidence  interval  obtained  when  calculating  the  mean  direction  from  the 
combined  DD-CM  and  GD-VM  results,  plus  the  mean  of  the  standard-deviation  ellipses  in 
Figure  4.20a,  plus  the  mean  of  the  standard-deviation  ellipses  in  Figure  4.20b  (calculated  using 
the  Fisher  precision  parameters  k  for  each).  Error  ellipses  are  not  calculated  for  the  two 
different  minimum  error  solutions  for  area  2  from  the  DD-CM  algorithm. 
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All  northern  paleopoles  reversed  into  southern  hemisphere 


A  SPA  negative  gravity  anomalies  A  Topography  paleopole 

(Garrick-Bethell  et  al.  2014) 


(C)  (D) 


Arkani-Hamed  et  al.  (2014),  isolated  anomalies  Arkani-Hamed  et  al.  (2014),  demagnetized  crust 
Figure  4.22.  Paleopoles  from  SPA  anomalies  compared  with  other  published  paleopoles  from 
other  anomalies.  Both  figures  show  the  southern  hemisphere  (all  points  are  southern  latitudes). 
(A)  Our  results  from  Figure  4.20,  with  all  northern  poles  reversed  into  the  southern  hemisphere. 
The  mean  paleopole  obtained  from  all  data  at  area  2  (see  Figure  4.21)  has  been  omitted  in  this 
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figure.  Red  and  blue  triangles  are  possible  paleopoles  for  SPA  immediately  after  its  formation, 
assuming  an  extreme  range  of  density  anomaly  models;  red  (blue)  points  are  positive  (negative) 
density  anomalies  (see  Section  4.5. 3 A).  SPA’s  center  (a  limiting  paleopole  for  very  negative 
density  models  of  SPA)  is  at  (-53.2°  S,  191°  E)  (black  triangle).  The  magenta  triangle  shows  the 
paleopole  dervied  from  the  tidal  component  of  the  Moon’s  topography,  outside  of  large  basins 
[Garrick-Bethell  et  al,  2014a].  (B)  Paleopoles  from  the  11  anomalies  reported  by  [Takahashi  et  al, 
2014].  No  anomalies  are  within  SPA.  Blue  (red)  points  represent  inversions  from  Kaguya  (Lunar 
Prospector)  data.  In  most  cases,  multiple  points  represent  paleopoles  from  single  sites.  Stars 
represent  the  means  of  the  two  clusters.  (C)  Paleopoles  from  10  isolated  anomalies  reported  by 
Arkani-Hamed  and  Boutin  (2014).  (D)  Paleopoles  from  the  crust  at  20  impact  craters  reported 
by  Arkani-Hamed  and  Boutin  (2014). 
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Figure  4.23.  Area  1  best  fit  results,  day  123,  2009.  Figure  details  as  in  Figure  4.9. 


0  < 

■-ioi 


117 


Lat  (deg  N)  Lat  (deg  N)  Lat  (deg  N)  Lat  (deg  N) 


LP  Data 


DD-CM 


GD-VM 


-23 

-24 

-25 

-26 

-27 


Area  1 ;  31 .1  km 


-186  -184  -182 


-186  -184  -182  -186  -184  -182 

Longitude  (deg  E)  Day  061 ,  1999 


Figure  4.24.  Area  1  best  fit  results,  day  61, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.25.  Area  1  best  fit  results,  day  142, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.26.  Area  2  best  fit  results,  day  69,  2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.27.  Area  2  best  fit  results,  day  96,  2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.28.  Area  2  best  fit  results,  day  33, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.29.  Area  3  best  fit  results,  day  69, 2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.30.  Area  3  best  fit  results,  day  96,  2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.31.  Area  3  best  fit  results,  day  61, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.32.  Area  4  best  fit  results,  day  151,  2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.33.  Area  4  best  fit  results,  day  115, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.35.  Area  5  best  fit  results,  day  96, 2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.36.  Area  5  best  fit  results,  day  151,  2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.37.  Area  5  best  fit  results,  day  142, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.38.  Area  6  best  fit  results,  day  96,  2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.39.  Area  6  best  fit  results,  day  33, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.40.  Area  6  best  fit  results,  day  115, 1999.  Figure  details  as  in  Figure  4.9. 


LP  Data 


cn 

0 

■D 


-22 


-24 


_  H 


DD-CM 


GD-VM 


20  \ 
15  " 
10  ; 
5  t 


U) 

0 

■D 


-22 


-24 


cn 

0 

■D 


-22 


-24 


• 

• 

1 

'*1  ^ 

•  • 

y 

10 


i-10 


-168  -166  -164  -162 

Area  6;  22.6  km 


-168  -166  -164  -162 

Longitude  (deg  E) 


-168  -166  -164  -162 

Day  142,  1999 


D 


10 

0 

-10 


Figure  4.41.  Area  6  best  fit  results,  day  142, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.42.  Area  7  best  fit  results,  day  96,  2009.  Figure  details  as  in  Figure  4.9. 


Figure  4.43.  Area  7  best  fit  results,  day  33, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.44.  Area  7  best  fit  results,  day  115, 1999.  Figure  details  as  in  Figure  4.9. 


Figure  4.45.  Area  8  best  fit  results,  day  96,  2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.46.  Area  8  best  fit  results,  day  151,  2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.47.  Area  8  best  fit  results,  day  61, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.48.  Area  8  best  fit  results,  day  142, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.49.  Area  9  best  fit  results,  day  96,  2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.50.  Area  9  best  fit  results,  day  61, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.51.  Area  9  best  fit  results,  day  88, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.52.  Area  9  best  fit  results,  day  115, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.53.  Area  9  best  fit  results,  day  142, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.54.  Area  10  best  fit  results,  day  151,  2009.  Figure  details  as  in  Figure  4.9. 
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Figure  4.55.  Area  10  best  fit  results,  day  61, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.56.  Area  10  best  fit  results,  day  88, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.57.  Area  10  best  fit  results,  day  142, 1999.  Figure  details  as  in  Figure  4.9. 
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Figure  4.58.  For  Area  1,  source  model  obtained  with  GD-VM  algorithm  [cf.  Hemingway  and 
Garrick-Bethell,  2012].  Each  square  represents  a  single  dipoles  covering  0.25  x  0.25  degrees 
(~5xl0’ m^);  colors  indicate  each  dipole’s  total  magnetic  moment.  Number  of  gridded  dipoles  is 
shown  in  Table  4.1.  Dataset  is  the  same  as  in  Figure  4.8. 
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Figure  4.59.  For  Area  2,  source  model  obtained  with  GD-VM  algorithm.  Details  are  as  in  Figure 
4.58,  except  that  the  dataset  is  the  same  as  in  Figure  4.9. 
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Figure  4.60.  For  Area  3,  source  model  obtained  with  GD-VM  algorithm.  Details  are  as  in  Figure 
4.58,  except  that  the  dataset  is  the  same  as  in  Figure  4.10. 
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Figure  4.61.  For  Area  4,  source  model  obtained  with  GD-VM  algorithm.  Details  are  as  in  Figure 
4.58,  except  that  the  dataset  is  the  same  as  in  Figure  4.11. 


136 


Longitude  (degrees  E) 

Figure  4.62.  For  Area  5,  source  model  obtained  with  GD-VM  algorithm.  Details  are  as  in  Figure 
4.58,  except  that  the  dataset  is  the  same  as  in  Figure  4.12. 


Figure  4.63.  For  Area  6,  source  model  obtained  with  GD-VM  algorithm.  Details  are  as  in  Figure 
4.58,  except  that  the  dataset  is  the  same  as  in  Figure  4.13. 
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Figure  4.64.  For  Area  7,  source  model  obtained  with  GD-VM  algorithm.  Details  are  as  in  Figure 
4.58,  except  that  the  dataset  is  the  same  as  in  Figure  4.14. 
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Figure  4.65.  For  Area  8,  source  model  obtained  with  GD-VM  algorithm.  Details  are  as  in  Figure 
4.58,  except  that  the  dataset  is  the  same  as  in  Figure  4.15. 
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Figure  4.66.  For  Area  9,  source  model  obtained  with  GD-VM  algorithm.  Details  are  as  in  Figure 
4.58,  except  that  the  dataset  is  the  same  as  in  Figure  4.16. 
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Figure  4.67.  For  Area  10,  source  model  obtained  with  GD-VM  algorithm.  Details  are  as  in  Figure 
4.58,  except  that  the  dataset  is  the  same  as  in  Figure  4.17. 
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Chapter  5 


Phase-Dependent  Atmospheric 
Retrievals  on  Gas  Giant  Planets  in 
Reflected  Light 


This  chapter  is  a  modified  reprint  of  M.  Nayak,  R.  Lupu,  M.  Marley,  J.  Fortney,  T. 
Robinson  and  N.  Lewis  (2016),  Atmospheric  Retrieval  for  Direct  Imaging 
Spectroscopy  of  Gas  Giants  In  Reflected  Light  II:  Orbital  Phase  and  Planetary 
Radius,  accepted.  Publications  of  the  Astronomical  Society  of  the  Pacific. 

5.1  Abstract 

Future  space-based  telescopes,  such  as  the  Wide-Field  Infrared  Survey 
Telescope  (WFIRST),  will  observe  the  reflected-light  spectra  of  directly  imaged 
extrasolar  planets.  Interpretation  of  such  data  presents  a  number  of  novel  challenges, 
including  accounting  for  unknown  planet  radius  and  uncertain  stellar  illumination 
phase  angle.  Here  we  report  on  our  continued  development  of  Markov  Chain  Monte 
Carlo  retrieval  methods  for  addressing  these  issues  in  the  interpretation  of  such  data. 
Specifically  we  explore  how  the  unknown  planet  radius  and  potentially  poorly  known 
observer-planet-star  phase  angle  impacts  retrievals  of  parameters  of  interest  such  as 
atmospheric  methane  abundance,  cloud  properties  and  surface  gravity.  As  expected. 
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the  uncertainty  in  retrieved  values  is  a  strong  function  of  signal-to-noise  ratio  (SNR) 
of  the  observed  spectra,  particularly  for  low  metallicity  atmospheres,  which  lack  deep 
absorption  signatures.  Meaningful  results  may  only  be  possible  above  certain  SNR 
thresholds;  for  cases  across  a  metallicity  range  of  1-50  times  solar,  we  find  that  only 
an  SNR  of  20  systematically  reproduces  close  to  the  correct  methane  abundance  at  all 
phase  angles.  However,  even  in  cases  where  the  phase  angle  is  poorly  known  we  find 
that  the  planet  radius  can  be  constrained  to  within  a  factor  of  two.  We  find  that 
uncertainty  in  planet  radius  decreases  at  phase  angles  past  quadrature,  as  the  highly 
forward  scattering  nature  of  the  atmosphere  at  these  geometries  limits  the  possible 
volume  of  phase  space  that  relevant  parameters  can  occupy.  Finally,  we  present  an 
estimation  of  possible  improvement  that  can  result  from  combining  retrievals  against 
observations  at  multiple  phase  angles. 

5.2  Introduction 

Transit  and  radial  velocity  (RV)  surveys  have  been  highly  successful  in 
detecting  short-period  exoplanet  systems,  and  have  allowed  the  compilation  of  a 
statistical  picture  of  the  bulk  properties  of  inner  planetary  systems.  However,  the  next 
frontier  in  exoplanet  studies  is  space-based  direct  imaging  and  spectroscopy  using 
optical-wavelength  telescopes,  coronagraphs  and  integral  field  spectrographs.  Such 
instruments  will  allow  the  characterization  of  colder  or  self-luminous  planets  that 
orbit  farther  from  their  parent  star.  The  upcoming  Wide-Field  Infra-Red  Survey 
Telescope  (WFIRST)  space  telescope  will  feature  a  space-based  high-contrast 
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coronagraph  for  imaging  and  spectroscopic  studies  of  planets  around  nearby  stars 
(Spergel  et  al.  2013).  It  will  perform  spectroscopy  of  extrasolar  planets  in  reflected 
light  at  spectral  resolutions  (R)  of  R~70,  in  wavelengths  ranging  from  -600-970  nm. 
Unlike  transit  spectroscopy,  which  is  able  to  probe  to  atmospheric  pressures  of  -1 
mbar  (Kreidberg  et  al.  2014),  or  ~1  bar  in  combination  with  emission  spectra,  direct 
imaging  has  the  potential  to  probe  deeper  into  the  atmosphere,  to  the  pressure  at 
which  atmospheric  aerosols  become  optically  thick  [Morley  et  al,  2015]. 

To  support  the  definition  of  future  direct  imaging  missions  and  to  enhance 
their  science  returns,  we  have  been  developing  a  set  of  tools  to  characterize  gas  giant 
planet  atmospheric  and  physical  properties  using  reflected  light  spectroscopy,  given 
anticipated  instrument  parameters  from  WFIRST.  In  Lupu  et  al.  (2016),  the  first  in 
this  series,  retrievals  of  atmospheric  methane  abundances  and  basic  cloud  properties 
using  Markov  Chain  Monte  Carlo  (MCMC)  techniques  were  explored,  assuming 
planets  with  known  radii  were  observed  at  full  phase.  [Lupu  et  al,  2016]  built  on 
previous  work  by  members  of  our  group  to  create  models  of  reflected  light  spectra, 
beginning  with  Marley  et  al.  (1999)  and  leveraging  albedo  variations  as  a  function  of 
cloud  structure,  mass,  metallicity,  planet  phase  and  star-planet  separation  by  Cahoy  et 
al.  (2010).  Other  contributions  in  this  field  have  included  Sudarsky  et  al.  (2000,  2003) 
and  Burrows  et  al.  (2004).  All  these  studies  of  reflected  light  spectra  of  exoplanets 
modeled  the  planets  at  full  phase  [Marley  et  al,  1999;  Lupu  et  al,  2016],  an  average 
phase  [Sudarsky  et  al,  2003]  or  at  a  set  of  specified  phase  angles  [Sudarsky  et  al. 


142 


2005;  Cahoy  et  al,  2010]  and  implicitly  assumed  that  the  incident  flux  and  planet  size 
were  known. 

However  during  a  real  observation  campaign,  several  factors  that  control  the 
reflected  flux  will  be  poorly  known.  First,  planets  will  be  observed  at  a  variety  of 
different  points  along  their  orbits.  Depending  on  the  fidelity  with  which  the  orbit  is 
constrained,  the  star-planet-observer  angle  (phase  angle),  and  thus  the  illumination 
phase,  may  not  be  well  known.  The  instantaneous  distance  of  the  planet  from  its  star, 
and  thus  the  incident  flux,  will  almost  certainly  not  be  perfectly  known.  Likewise, 
planet  radii  will  not  be  constrained,  except  by  the  observed  brightness  of  the  planet  at 
a  variety  of  wavelengths  and  the  mass-radius  relationship  for  those  planets  with 
masses  constrained  by  radial  velocity  measurements.  Any  uncertainties  in  orbital 
phase  will  further  obscure  planet  radius  determination. 

Figure  5.1  illustrates  the  degenerate  nature  of  planet  radius  with  increasing 
planet  phase  in  scattered  light;  the  brightness  of  a  planet  can  decrease  either  with 
decreasing  planet  radius  or  increasing  phase.  In  other  words,  a  large  planet  at 
quadrature  (phase  angle  a  =  90°)  and  a  smaller  planet  at  full-phase  (a  =  0°)  could  not 
be  distinguishable  solely  by  their  relative  brightness.  As  we  shall  show,  this 
degeneracy  has  a  significant  effect  on  the  quality  of  the  resulting  retrievals  on  other 
parameters  of  interest,  including  methane  abundance  and  cloud  properties.  In  this 
work  we  explore  the  bounds  of  these  mutually  dependent  parameters  and  determine 
signal-to-noise  (SNR)  requirements  for  scientifically  interesting  conclusions  to  be 
drawn. 
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The  goal  of  the  current  work  is  therefore  to  characterize  the  effects  of 
changing  planet  phase  on  retrievals  of  atmosphere  molecular  abundances,  surface 
gravity  and  cloud  properties  when  planetary  radius  is  unknown.  The  paper  is 
organized  as  follows:  Section  5.3  provides  background  on  our  approach  to  modeling 
the  phase  angle  and  clouds,  our  MCMC  formulation  and  our  chosen  exoplanet  test 
cases  of  HD  99492c  and  HD  192310c;  Section  5.4  details  MCMC  retrieval  results 
and  the  use  of  posterior  probability  plots  to  extract  68%  confidence  intervals  on 
parameters  of  interest.  Section  5.5  contains  our  discussion  of  how  planet  radius,  phase 
angle,  methane  abundance  and  cloud  properties  were  constrained  in  the  presence  of 
planet  phase  and  radius  uncertainties;  finally.  Section  5.6  presents  our  conclusions. 

5.3  Background 

In  this  section,  we  provide  a  brief  overview  of  some  key  concepts  discussed  in 
the  paper.  This  work  builds  on  previous  work  by  several  authors  to  create  models  of 
albedo  spectra  of  extrasolar  giant  planets;  a  thorough  description  may  be  found  in 
[Lupu  et  al,  2016].  Our  initial  study  reported  in  that  paper  represented  the  first  time 
molecular  abundances  and  cloud  properties  were  simultaneously  retrieved  using 
Bayesian  inference  tools  applied  to  simulated  scattered  light  spectra  of  cool  extrasolar 
giant  planets.  Similar  applications  of  Bayesian  methods  to  exoplanet  studies  include 
work  by  Irwin  et  al.  (2008)  and  Benneke  &  Seager  (2012),  among  others. 

5.3.1  Albedo  Model 
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To  compute  the  thermal  structure  of  each  model  planet’s  atmosphere,  we  use  a 
ID  radiative-convective  equilibrium  model  based  on  that  developed  for  Titan 
(McKay  et  al.  1989)  and  solar  system  giant  planets  and  exoplanets  [Marley  and 
McKay,  1999;  Marley  et  al,  1999].  The  methane  opacity  at  optical  wavelengths  is 
taken  from  Karkoschka  (1994)  and  collision-induced  absorption  from  Freedman  et  al. 
(2008).  In  this  paper  our  test  planets  are  cold  and  we  neglect  H2O  opacity.  Given  a 
self-consistent  model  atmospheric  profile,  we  compute  an  albedo  spectrum  following 
the  methods  described  in  Cahoy  et  al.  (2010).  Cloud  scattering  is  treated  with  a  two- 
term  Heyney-Greenstein  function,  which  captures  moderate  backscattering  and  high 
forward  scattering  as: 

Vtt-hg  “  Pncid’  +  ^Phg  (~f  <  ^)  ^  ^ 

where: 


PHG(g>G)  = 


(l+g^-2gcos6)^-^ 


(5.2) 


Here,  g  is  the  scattering  asymmetry  factor,  a  measure  of  the  preferred 
direction  of  light  scattered  by  aerosol  particles;  this  is  a  retrievable  quantity. 
Integrating  over  the  emergent  intensity  using  Gaussian-Chebyshev  quadrature  {[Lupu 
et  al.,  2016],  originally  from  Horak  (1950)  and  Horak  &  Little  (1965))  yields  model 
albedo  spectra  for  the  planet.  In  this  paper  we  treat  the  phase  angle  as  fully  unknown. 
We  allow  it  to  vary  from  full  phase  (a  =  0°)  to  a  =  135°,  at  which  angle  the  flux  for  a 
Lambertian  sphere  is  <5%  of  that  at  full  phase.  In  reality,  we  will  have  some 
constraint  on  phase  angle  from  radial  velocity  (RV)  observations,  which  will  improve 
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retrievals;  our  results  may  therefore  be  treated  as  a  worst-case  scenario  (but  see 


Section  4.5.2).  Similarly,  even  though  observations  at  full  phase  are  not  possible  for 
direct  imaging,  we  include  that  possibility  to  understand  how  the  quality  of  retrievals 
at  other  phases  may  relate  to  those  at  zero  phase,  which  were  elucidated  in  [Lupu  et 
al,  2016]. 


5.3.2  MCMC  Methodology 

For  our  implementation  of  Markov  Chain  Monte  Carlo  (MCMC)  methods,  we 
follow  the  approach  developed  for  massively  parallel  implementations  by  Lupu  et  al. 
(2016)  and  use  emcee,  an  open-source  affine  invariant  ensemble  MCMC  sampler 
(Foreman-Mackey  et  al.  2013;  Goodman  &  Weare  2010).  Given  a  set  of  well-chosen 
bounds  on  retrievable  parameters  (“priors”),  this  approach  efficiently  samples  the 
parameter  space  and  allows  for  massively  parallel  computation.  For  each  retrievable 
parameter,  this  implementation  employs  multiple  MCMC  chains  in  parallel.  Full 
details,  and  a  comparison  of  emcee  to  another  multimodal  nested  sampling  algorithm 
(MultiNest),  are  discussed  in  [Lupu  et  al.,  2016]. 

Table  5.1.  Description  of  retrievable  parameters  in  two-layer  cloud  model  [Marley  et  al.,  2014],  as 

well  as  priors  used  for  MCMC  runs. 


Quantity 

Description 

Priors 

Descriptor 

fcH4 

Molecular  abundance  of  Methane 

-8  to  -2 

log  space 

g 

Surface  acceleration  due  to  gravity 

1-300 

m/s^ 

R 

Planet  Radius 

0.1-100 

Jupiter  radius 

dPi 

Pressure  difference:  Top  of  lower  cloud  to  Bottom  of  upper  cloud 

-2  to  2 

log,  bar 

dP, 

Pressure  difference:  Bottom  of  upper  cloud  to  Top  of  upper  cloud 

-2  to  2 

log,  bar 

T 

Total  optical  depth,  upper  cloud 

-2  to  2 

log,  bar 

nxi 

Single  scattering  albedo,  upper  cloud 

0.01  to  0.9999 

g 

Asymmetry  Factor 

0.01  to  0.9999 

P 

Pressure,  top  of  lower  cloud 

-2  to  1.5 

log,  bar 

TD-2 

Single  scattering  albedo,  lower  cloud 

0.01  to  0.9999 

0 

Planet  phase  angle 

0  to  135 

degrees 
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As  also  described  in  {Lupu  etal,  2016],  we  apply  these  methods  to  simulated 
spectral  datasets,  computed  as  in  Section  4.3.1,  to  retrieve  quantities  of  interest  in  the 
interpretation  of  exoplanet  spectra.  Table  5.1  lists  all  eleven  retrievable  quantities  that 
are  estimated  by  our  Markov  Chain  Monte  Carlo  (MCMC)  routine.  It  assumes  that 
the  atmosphere’s  major  absorber  is  solely  methane,  with  H2-He  background  gas. 
Priors  are  set  to  allow  values  to  range  across  six  orders  of  magnitude  for  methane 
abundance,  2.5  orders  of  magnitude  for  surface  gravity  and  three  orders  of  magnitude 
for  planetary  radius  (Table  5.1).  Of  course,  these  are  extremely  large  ranges;  in  reality 
better  constraints  are  expected.  For  example,  astrometry  combined  with  RV 
constraints  will  likely  determine  the  planet  mass  to  within  a  factor  of  two.  Likewise, 
the  mass-radius  relationship  trivially  demonstrates  that  a  Jupiter  mass  planet  would 
never  have  a  radius  of  100  Rj.  However,  the  exploration  of  a  large  parameter  space 
can  be  valuable  in  permitting  a  greater  number  of  feasible  solutions,  enabling  a  better 
understanding  of  the  degeneracies  inherent  in  the  problem.  As  an  example,  all  else 
being  equal,  lower  methane  abundance  at  lower  gravity  can  produce  a  similar 
absorption  feature  to  a  higher  methane  abundance  at  a  higher  surface  gravity  [Marley 
et  al,  2014;  Lupu  et  al,  2016];  the  use  of  MCMC  to  explore  this  large  parameter 
space  may  be  useful  to  probabilistically  distinguish  between  the  two  cases. 

Also  among  the  retrievable  parameters,  exoplanet  cloud  and  haze  aerosols  are 
parameterized  using  an  improved  version  of  the  simple  two-layer  cloud  model  first 
detailed  in  Marley  et  al.  (2014).  Using  the  two-layer  cloud  model,  illustrated  in  Figure 
5.2,  we  create  a  model  cloud  and  a  noise-free  albedo  spectrum  (Marley  et  al.  2014). 
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These  quantities  are  also  used  in  [Lupu  et  al,  2016],  which  contrasts  the  two-layer 
model  against  simpler  one-layer  and  no  cloud  models. 

To  understand  how  instrumental  and  astrophysical  parameters  can  affect 
observed  spectra,  and  include  these  effects  in  our  retrievals,  we  apply  a  noise  model 
to  the  noise-free  spectrum,  which  includes  convolution  with  an  instrument  point 
spread  function  (PSF)  to  an  appropriate  spectral  resolution,  notionally  representative 
of  that  expected  for  WFIRST.  Parameters  of  the  noise  model  are  presented  in  Table 
5.2  and  the  implementation  follows  Robinson  et  al.  (2015).  Each  simulated  data  point 
is  drawn  from  a  normal  distribution,  with  the  mean  given  by  the  planet-star  flux  ratio 
and  the  standard  distribution  given  by  the  noise  model.  This  noisy  spectrum  then 
becomes  the  input  to  the  MCMC  retrieval  code. 

We  make  three  notes  here:  firstly,  at  the  time  of  this  writing,  the  WFIRST 
coronagraph  instrument  parameters  are  still  under  active  study  and  refinement 
[Harding  et  al,  2015].  However,  our  adopted  values  (Table  5.2)  are  meant  to  be 
representative  of  the  noise  levels  that  the  mission  is  expected  capable  of  achieving,  as 
it  is  understood  to  be  in  mid-2016.  Secondly,  the  instrument-representative  approach 
employed  here  differs  from  [Lupu  et  al,  2016],  which  uses  a  general  synthetic  noise 
model.  Finally,  here  we  do  not  attempt  to  retrieve  the  atmospheric  temperature- 
pressure  profile,  as  the  reflected  light  spectra  are  only  weakly  dependent  on  the 
profile.  Variations  in  gravity  do  alter  the  scale  height  and  atmospheric  density  and 
these  efafects  are  accounted  for.  A  future  paper  in  this  series  will  explore  atmospheric 
temperature  profile  retrievals. 
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Table  5.2.  Parameters  used  in  the  notional  WFIRST  noise  model.  Details  on  implementation 

follow  [Robinson  etal.,  2015] 


Item 

WFIRST 

representative  value 

Unit 

Dark  current 

5.00E-04 

s-1 

Telescope  diameter 

2.4 

m 

Read  noise 

0.2 

per  pixel 

System  throughput 

0.037 

Angular  size  of  lenslet 

0.017 

arcsecond 

Inner  working  angle 

2.7 

X/D 

Outer  working  angle 

10 

X/D 

Size  of  photometric  aperture 

1.5 

X/D 

Contrast  floor 

l.OOE-10 

Minimum  wavelength 

0.6 

urn 

Maximum  wavelength 

0.95 

urn 

Spectral  resolution 

R  =  70 

Here  we  do  not  attempt  to  retrieve  the  atmospheric  temperature-pressure 
profile,  as  the  reflected  light  spectra  are  only  weakly  dependent  on  the  profile. 
Variations  in  gravity  do  alter  the  scale  height  and  atmospheric  density  and  these 
effects  are  accounted  for.  Future  work  will  explore  atmospheric  profile  retrievals. 


5.3.3  Test  Cases:  Synthetic  HD  99492c  (Planet  A)  and  HD 
192310c 

Four  test  cases  across  two  idealized  exoplanets  encompass  our  efforts  to 
contrast  the  relative  effects  of  changing  planet  mass.  Represented  in  the  retrievable 
parameters  by  surface  gravity  and  planet  radius,  the  planet  mass  in  turn  controls  other 
retrieved  properties.  For  a  non-solar  system  test  case,  [Lupu  et  al,  2016]  used  HD 
99492c  [Marcy  et  al,  2005];  our  first  test  case  is  also  inspired  by  HD  99492c.  The 
inferred  mass  of  this  planet  is  0.36+0.02  Mj  (Meschiari  et  al.  2010,  Table  2);  at  a 
semi-major  axis  of  5.4  AU  from  its  star,  models  from  Fortney  et  al.  (2007)  suggest  a 
radius  for  this  planet  of  -0.9  Rj.  Kane  et  al.  (2016)  report  that  HD  99492c  is  in  fact  an 
artifact  attributable  to  variability  of  the  host  star  and  not  a  planet.  Therefore  we  treat 
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spectra  generated  for  HD  99492c  as  a  synthetic  case  study  to  baseline  our  results 
against,  for  a  planet  almost  of  Jupiter  radius  with  a  methane-dominated  atmosphere. 
To  permit  observations  at  all  phase  angles  to  be  above  the  WFIRST  contrast  floor 
(Table  5.2),  we  “relocate”  our  synthetic  planet  to  2  AU  from  its  star  and  5  parsecs 
from  the  telescope;  the  actual  values  for  HD  99492c  are  5.4  AU  /  18  parsecs  [Marcy 
et  al,  2005].  We  refer  to  this  synthetic  HD  99492c  analog  as  “Planet  A”  for  the 
remainder  of  this  work. 

The  second  test  case  explores  a  planet  with  an  order  of  magnitude  smaller 
mass  (Figure  5.3).  HD  192310c,  also  known  as  Gliese  785  c,  has  M  sin  i  of 
0.076+0.016  Mj  [Pepe  et  al,  2011].  Using  mass-radius  relationships  from  Fortney  et 
al.  (2007)  we  infer  a  radius  of  0.75  Rj.  Using  the  notional  parameters  in  Table  5.2, 
this  planet  will  be  near  the  coronagraph  inner  working  angle  (IWA)  at  600 
nm.  However,  since  the  WFIRST  coronagraphs  are  still  under  development,  their  in¬ 
flight  performance  may  yet  change.  Thus,  we  choose  not  to  relocate  the  planet  as  for 
HD  99492c  /  Planet  A,  and  instead  use  HD  192310c  as  a  test  case  for  situations  where 
a  planet  lies  near  enough  to  the  coronagraph  IWA  to  cause  the  noise  to  be  dominated 
by  stellar  leakage.  Using  this  approach,  we  produce  three  distinct  models  of  HD 
192310c  at  metallicities  of  lx,  lOx  and  50x  times  solar  abundance,  which  produce 
successively  deeper  methane  features  (see  Figure  5.6).  Figure  5.4  contrasts  model 
spectra  for  Planet  A  and  the  50x  solar  metallicity  HD  192310c  at  a  =  0°  and  a  =  90°. 
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5.4  Results 

5.4.1  Retrieved  best-fit  spectra  by  MCMC 


We  generate  model  albedo  spectra  for  Planet  A  (Figure  5.5)  and  HD  192310c 
(Figure  5.6)  at  varying  phase  angles  similar  to  Figure  5.4;  in  this  study  we  explore 
phase  angles  of  30°,  40°,  60°,  70°,  90°  and  120°.  Both  resulting  spectra  are  then 
combined  with  an  instrument-specific  model,  in  this  case,  the  notional  WFIRST  noise 
model  in  Table  5.2  [Robinson  et  ah,  2015]. 

While  [Lupu  et  al,  2016]  performs  retrievals  against  an  albedo  spectrum, 
since  this  work  is  also  concerned  with  planet  radius,  we  retrieve  against  the  planet-to- 
star  flux  ratio,  or  contrast  spectrum.  Contrast  spectra  are  created  for  a  range  of  signal- 
to-noise  ratios  (SNR)  to  explore  observational  limits  and  their  effect  on  observations. 
Here  we  define  the  SNR  to  be  centered  at  A=  0.6  micron  with  an  8.6  nm  wide 
bandpass.  SNRs  of  5,  10  and  20  are  explored,  as  in  [Lupu  et  al,  2016].  To  be  clear, 
since  SNR  is  wavelength-dependent,  these  values  refer  to  the  SNR  at  0.6  pm. 

For  each  of  the  eleven  retrievable  parameters  (Table  5.1),  we  employ  twelve 
MCMC  chains,  or  “walkers”,  for  a  total  of  132  chains.  Each  chain  was  then  run  for 
2500  iterations  for  a  final  sample  chain  of  330,000  samples.  Selected  cases  were  run 
for  4000  iterations  to  ensure  that  the  MCMC  algorithm  did  not  get  stuck  in  a  local 
minimum  and  was  exploring  the  entire  parameter  space;  returned  ranges  were  found 
to  be  nearly  identical  to  the  2500  iteration  run,  so  we  restrict  ourselves  to  2500 
iterations  for  all  the  results  presented  here.  We  identified  the  best  fit,  la  and  2a-range 
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of  retrieved  spectra  for  Planet  A;  the  retrieved  models  match  the  “true”  spectra  well. 
Similar  excellent  fits  are  seen  for  phase  angles  between  30°  and  120°  despite 
decreasing  contrast  signals  at  larger  phase  angles  (Figure  5.7). 

We  perform  a  similar  study  against  HD  192310c  using  the  same  MCMC 
parameters  as  for  Planet  A.  For  this  planet,  we  generate  three  separate  test  cases  by 
constructing  forward  models  and  performing  retrievals  at  three  different  metallicities, 
namely,  one,  ten  and  fifty  times  solar  values  (lx,  lOx,  50x).  For  each  metallicity  case, 
we  generate  model  albedo  spectra  at  30°,  40°,  60°,  70°,  90°  and  120°  phase  angle  and 
apply  the  notional  WFIRST  noise  model  to  them. 

The  resulting  spectra  reveal  a  variety  of  methane  absorption  signatures  (Figure 
5.6).  Notably,  because  of  the  relatively  high  cloud,  the  lx  solar  case  exhibits 
particularly  subdued  methane  absorption  features.  Figure  5.8  -  Figure  5.10  illustrate 
the  corresponding  spectral  recoveries  across  three  values  of  signal-to-noise  (SNR) 
and  three  values  of  metallicity,  at  phase  angles  of  30°  and  90°. 

5.4.2  Inferring  Phase-Dependent  Relationships  from 
Posterior  Probability  Distributions 

We  assemble  retrieval  results  similar  to  those  shown  in  Section  5.4.1  into 
posterior  probability  plots,  which  graphically  show  the  marginal  probability 
distribution  between  every  retrieved  parameter  pair.  Figure  5.11  is  an  example  of  this 
plot.  Here,  darker  colors  represent  a  higher  probability  of  the  solution  lying  in  that 


152 


region,  and  the  diagonal  of  the  plot  shows  the  marginalized  probability  distribution,  to 
a  68%  confidence  interval,  for  each  retrievable  parameter  in  Table  5.1. 

A  broad  distribution  means  that  the  parameter  is  largely  unconstrained,  as  is 
the  case  for  phase  angle  (Figure  5.11,  marker  a)  or  surface  gravity.  Conversely,  a 
sharp  peak  in  the  distribution  and  small  ranges  on  returned  values  means  that  the 
parameter  can  be  well  determined,  such  as  planet  radius  (Figure  5.11,  marker  b).  In 
other  cases,  more  general  relationships  can  be  inferred,  for  example,  lower  limits  to 
the  albedo  of  the  top  cloud  (Figure  5.11,  marker  c)  and  methane  abundance.  Such 
plots  illuminate  how  variations  in  retrieved  values  vary  with  changes  in  the  "true" 
planet  phase,  as  well  as  interrelationships  between  other  parameters.  A  brief 
discussion  of  general  trends  apparent  in  the  relationships  of  phase  angle  with  other 
retrievable  parameters  follows;  though  we  use  the  Planet  A  case  (Figure  5.11)  to 
highlight  these  trends,  they  are  also  seen  for  HD  192310c. 

Figure  5.12  focuses  on  selected  panels  from  Figure  5.1 1.  First,  the  relationship 
between  radius  and  phase  is  in  line  with  our  conceptual  understanding  from  Figure 
5.1:  for  larger  phase  angles  (i.e.  a  more  crescent  phase),  a  larger  planet  radius  is 
favored.  Even  at  a  relatively  high  phase  angle  (60°)  a  clear  detection  of  methane,  with 
a  lower  limit  to  the  atmospheric  mixing  ratio  of  larger  than  10'^  ^,  is  seen.  Surface 
gravity  appears  essentially  unconstrained,  although  as  we  will  show  later,  a  larger 
SNR  or  smaller  phase  angle  does  narrow  the  probable  range.  The  difficulty  in 
deriving  meaningful  constraints  on  gravity  from  reflection  spectra  is  discussed  in 
more  detail  in  \Lupu  et  al,  2016]. 
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The  MCMC  analysis  appears  to  constrain  the  top  cloud  well  (quantity  dP2,  see 
Table  5.1),  but  is  indeterminate  on  the  pressure  “gap”  between  the  cloud  layers 
(quantity  dPi).  The  pressure  of  the  bottom  cloud  (quantity  P2)  is  not  tightly 
constrained,  but  higher  probability  values  are  distributed  around  the  true  value.  This 
could  imply  either  that  the  bottom  cloud  is  not  well  constrained,  or  perhaps  that  a 
one-layer  cloud  model  is  more  suitable  for  this  planet.  Previous  MCMC  simulations 
on  HD  99492c  reach  a  similar  conclusion  [Marley  et  al,  2014]. 

We  generate  marginalized  probability  distributions  at  seven  phase  angles  (0°  - 
120°)  and  three  SNR  values  (5,  10,  20).  For  the  case  of  SNR  =  20,  Figure  5.13  shows 
the  relationship  between  retrieved  phase  angle  and  planet  radius  for  phase  angles 
from  30°  to  120°.  While  the  probable  ranges  on  retrieved  planet  phase  angle  are  large, 
at  both  low  and  high  phase  angles,  the  MCMC  algorithm  retrieves  best-fit  values 
close  to  the  true  value.  However,  at  phase  angles  between  45°  and  90°,  the  probable 
phase  angle  values  stretch  across  most  of  the  phase  angle  solution  space,  making  it 
more  difficult  to  obtain  good  values.  Observations  with  comparable  SNR  and 
multiple  phase  angles  can  therefore  be  extremely  valuable  in  determining  both  orbital 
characteristics,  if  unknown  or  uncertain,  and  narrowing  down  planet  radius. 

Finally,  we  collate  summary  plots  of  the  retrieved  parameters  at  all  seven 
phase  angles  (0°  -  120°)  and  SNR  values  (5,  10,  20)  for  all  four  test  cases.  The 
retrieved  values  of  methane  abundance,  surface  gravity,  planet  radius,  recovered 
phase  angle  and  cloud  pressures  are  respectively  plotted  against  changing  SNR  and 
phase  angle  for  Planet  A  (Figure  5.14),  HD  192310c  lx  (Figure  5.15),  lOx  (Figure 
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5.16)  and  50x  (Figure  5.17)  cases.  Colors  represent  the  size  of  the  68%  confidence 
interval  values,  seen  for  each  parameter  on  the  diagonal  of  probability  distribution 
plots  such  as  Figure  5.11. 

Several  trends  are  apparent  from  the  summary  figures.  Generally  speaking,  we 
find  that  retrievals  at  higher  SNR  ratios  (SNR  =  20)  place  correct  constraints  on  the 
atmospheric  methane  abundance,  to  within  an  order  of  magnitude.  The  low  SNR=5 
case  identifies  the  presence  of  methane,  but  the  abundance  is  highly  uncertain  (four 
orders  of  magnitude),  as  many  combinations  of  cloud  top  pressure,  phase  angle,  and 
gravity  are  able  to  adequately  reproduce  the  noisy  data. 

In  general,  the  high  bright  clouds  found  in  this  case  seem  to  lead  the  retrievals 
to  favor  cloud  tops  deeper  in  the  atmosphere  than  in  the  forward  model,  with  lower 
brightness  compensated  by  larger  planetary  radii.  This  is  seen  for  all  metallicities, 
particularly  at  low  phase  angles.  The  lx  metallicity  case  with  the  weakest  methane 
features  clearly  presents  a  particular  challenge,  even  at  SNR  of  10,  as  the  methane 
abundance  is  nearly  unconstrained.  For  all  three  metallicities  considered,  only  the 
case  with  an  SNR  of  20  systematically  reproduces  close  to  the  correct  methane 
abundance  at  all  phase  angles.  We  discuss  these  findings  further  in  the  next  section. 

5.5  Discussion 

In  this  section,  we  discuss  how  well  the  planet  radius,  phase  angle,  methane 
abundance  and  cloud  properties  were  constrained  in  the  presence  of  planet  phase  and 
radius  uncertainties.  We  focus  the  discussion  in  this  paper  on  the  newly  introduced 
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uncertain  phase  and  radius  determinations,  as  the  abundance  and  cloud  properties 
were  the  focus  of  [Lupu  et  al,  2016]. 

5.5.1  Methane  and  Radius  Retrieval 

For  the  cases  considered  here  we  assumed  that  phase  angle  was  almost 
completely  unconstrained.  While  this  is  a  situation  unlikely  to  be  encountered  for  a 
real  planet,  it  is  a  difficult  bounding  case  worthy  of  additional  study.  We  find  that 
phase  angle  is  generally  not  well  constrained  from  a  single  observed  spectrum,  since 
planet  radius  and  cloud  height  trade  against  phase  angle,  resulting  in  large 
uncertainties  in  all  parameters.  Generally  speaking,  with  no  prior  knowledge  of 
orbital  parameters  from  radial  velocity,  the  most  we  can  confidently  tell  about  the 
phase  angle  from  retrievals  is  whether  it  is  high  or  low  (above  or  below  -90°). 

One  might  expect  the  planet  radius  solution  space  to  be  similarly  large,  but 
this  is  not  the  case.  Despite  being  given  an  impossibly  large  range  of  0.1  -  100  Rj, 
even  at  a  low  SNR  of  5,  the  MCMC  routine  typically  returns  a  solution  within  a 
factor  of  two  of  the  true  value  (-1  Rj).  Regardless  of  true  phase  angle,  the  MCMC 
algorithm  must  match  the  observed  flux  from  the  planet.  At  more  crescent  phases,  the 
atmosphere  is  highly  forward  scattering  and  molecular  bands  are  weak.  Since  clouds 
are  relatively  less  important  at  such  scattering  angles,  this  drastically  reduces  the 
number  of  free  parameters.  Consequently,  we  find  the  most  accurate  radius  retrievals 
at  the  highest  phase  angles. 

5.5.2  The  impact  of  a  known  phase  angle 
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Given  the  difficulty  in  retrieving  the  true  phase  angle  from  a  single  spectrum, 
we  also  explored  the  impact  of  a  better-constrained  phase  angle,  as  might  be  expected 
during  an  observational  campaign.  To  study  this  case  we  chose  the  lOx  metallicity 
case  for  HD  192310c,  at  a  favorable  SNR  of  20.  Retrievals  were  performed  on  this 
case  at  multiple  phase  angles,  given  a  +70°  restriction  from  the  true  value,  on  the 
possible  values  of  the  angle. 

As  seen  in  Figure  5.18,  this  does  not  significantly  improve  retrievals  of 
gravity,  cloud  properties  or  the  methane  abundance.  The  only  noticeable  difference  is 
that  phase  angle  knowledge  helps  constrain  the  radius  of  the  planet  better,  improving 
the  radius  determination  by  a  factor  of  two.  Given  a  proxy  value  for  planet  mass  from 
radial  velocity  measurements  (M  sin  /),  by  improving  knowledge  of  the  planet  radius, 
prior  knowledge  of  the  phase  angle  best  helps  improve  the  estimate  of  surface  gravity 
(M/r^),  though  this  is  not  seen  directly  from  gravity  retrievals. 

5.5.3  Applying  an  intersection  criterion  to  multiple 
observations 

It  has  been  shown  that  planet  radius  retrievals,  for  example,  improve  with 
increasing  phase  angle,  whereas  retrievals  for  top  cloud  pressure  improve  with 
decreasing  phase  angle.  Bounds  on  quantities  such  as  methane  abundance  and  cloud 
properties  vary  significantly  with  SNR.  Simultaneous  retrievals  on  observations  taken 
at  multiple  phase  angles  would  therefore  likely  hold  promise  for  narrowing  the 
solution  space  of  best-fit  models.  While  we  did  not  perform  simultaneous  retrievals. 
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we  present  a  preliminary  investigation  into  their  utility  by  imposing  an  intersection 
criterion. 

The  intersection  criterion  is  defined  as  in  set  theory:  for  sets  A  and  B,  the 
intersection  of  sets  (A  fl  B)  defines  that  set  which  contains  only  elements  of  A  that 
also  belong  to  B.  For  two  observations  taken  at  different  phase  angles,  and  separate 
retrievals,  applying  an  intersection  criterion  means  only  those  solutions  that  appear  in 
both  retrievals  are  considered  valid.  Here  we  use  “observation”  to  refer  to  the 
integration  time  needed  to  produce  one  complete  -600-970  nm  spectrum,  at  a  given 
SNR.  The  intersection  criterion  is  illustrated  in  Figure  5.19.  While  clearly  an 
estimation,  simultaneous  retrievals  against  combined  phase-varying  datasets  are 
planned  as  future  work. 

The  underlying  idea  is  that  intersection  of  multiple  observations  at  varying 
phase  angles  may  determine  a  more  likely  range  for  parameters  of  interest.  We  begin 
by  determining  which  combination  of  phase  angles  will  be  likely  to  improve 
retrievals  the  most.  This  analysis  focuses  on  the  lowest  SNR  case,  as  the  uncertainty 
on  retrieved  results  with  one  set  of  observed  spectra  is  the  highest,  and  multiple 
observations  are  likely  to  have  the  most  impact.  The  HD  192310c  lOx  case  is  chosen 
again  here,  although  a  similar  analysis  may  be  conducted  for  any  planet  targeted  as 
part  of  an  observational  campaign. 

Figure  5.20  shows  the  improvement  in  68%  confidence  intervals  achieved  by 
applying  the  intersection  criterion  to  all  phase  angle  solutions,  for  HD  192310c  lOx, 
at  an  SNR  of  5.  Here,  “improvement”  denotes  the  difference  between  1)  confidence 
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intervals  obtained  using  single  observations  and  2)  intervals  obtained  by  applying  an 
intersection  criterion  to  multiple  observations  taken  at  differing  phase  angles.  A  value 
of  zero  for  improvement  represents  one  of  two  cases;  either  1 )  there  is  no  overlap 
between  retrieved  values  at  different  phases,  and  the  solutions  must  be  considered 
independently,  or  2j  returned  solutions  are  completely  identical  at  both  phase  angles. 
In  either  case,  considering  multiple  observations  does  not  represent  an  improvement 
over  a  single  observation.  Conversely,  peaks  represent  cases  where  confidence 
intervals  were  tightened  by  applying  an  intersection  criterion  to  multiple 
observations,  i.e.,  we  can  estimate  which  combination  of  phase  angles  may  be  most 
helpful  to  reduce  uncertainty  in  retrieved  quantities.  For  example,  for  methane,  if  a 
first  observation  is  taken  at  a=70°  (x-axis),  a  subsequent  observation  at  a=30° 
(yellow  line)  would  improve  the  68%  confidence  interval  by  approximately  1.4  orders 
of  magnitude  (y-axis). 

Figure  5.20  shows  that  a  steadily  increasing  improvement  in  estimates  of 
planet  radius  can  be  expected  with  two  observations  at  phase  angles  that  exceed  45°. 
For  example,  for  one  observation  at  70°  and  another  at  120°,  the  estimate  for  planet 
radius  can  be  improved  by  a  factor  of  two,  a  significant  improvement  when 
considering  that  the  best-fit  solution  from  a  single  observation  was  already  within  a 
factor  of  two  of  the  true  solution.  Similarly,  confidence  intervals  on  surface  gravity 
may  be  improved  by  as  much  as  55  -  60  ms'^  if  observations  are  gathered  at  0°  and 
70°  phase  angle,  though  gains  of  >30  ms'^  are  possible  with  other  combinations.  Both 
these  cases  were  also  illustrated  in  Figure  5.19. 
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Similarly,  the  uncertainty  in  methane  composition  from  retrievals  against  an 
observation  taken  at  near-quadrature  can  be  driven  down  by  almost  two  orders  of 
magnitude  if  combined  with  a  low  phase  angle  observation,  but  in  the  absence  of  such 
an  observation,  may  still  be  reduced  by  0.8  orders  with  an  observation  at  45°  phase 
angle.  The  intersection  criterion  presents  a  way  to  estimate  trends  in  phase-varying 
behavior;  future  simultaneous  retrievals  and  joint  probability  distributions  created 
from  multiple  observations  will  quantify  the  exact  improvement. 

Such  improvements  will  of  course  be  reduced  with  increasing  SNR  (smaller 
probability  bounds,  greater  overlap),  and  with  more  than  two  observations.  During  an 
actual  space-based  observational  campaign,  operational  or  other  constraints  may  limit 
the  ability  to  observe  a  planet  at  a  favorable  viewing  geometry.  Therefore,  we  now 
estimate  the  likelihood  of  improvement  in  confidence  intervals. 

We  randomly  choose  three  phase-varying  observations,  and  compare  the 
intersection  criterion  result  with  that  of  a  single  observation,  also  randomly  chosen. 
Figure  5.21  plots  the  improvement  with  the  intersection  criterion,  for  every  possible 
three-observation  combination  of  the  seven  phase  angles  studied  here,  for  both  Planet 
A  and  the  lOx  HD  192310c  case.  As  expected,  for  either  planet,  the  improvement 
generally  does  not  exceed  an  order  of  magnitude  for  the  SNR  =  20  case.  However  for 
SNR  =  5,  the  improvement  is  significant  in  all  cases,  even  for  largely  invariant 
parameters  such  as  gravity  (marker  d-e),  but  particularly  for  cloud  parameters  (marker 
f-g)  and  methane  abundance.  Up  to  3.5  orders  of  magnitude  in  improvement  for  the 
methane  abundance  is  seen  for  Planet  A  (marker  h),  depending  on  the  phase  angle 
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combinations;  recall  that  for  this  case,  we  were  unable  to  do  much  better  than 


determine  the  presence  of  methane. 

Similar  results  are  seen  for  four  and  more  observations.  Such  plots  allow  us  to 
build  an  idea  of  trends  for  improvement  in  uncertainty  estimates.  These  trends  will  be 
important  for  mission  planning  for  WFIRST;  these  also  present  a  starting  point  for 
estimating  science  return  for  realistic  mission  scenarios,  where  data  is  available  from 
multiple  observations  at  different  phase  angles  and  low  SNR. 

5.6  Conclusions 

We  have  studied  a  number  of  retrievals  on  simulated  phase-varying  spectra, 
incorporating  different  metallicities,  star-planet  fluxes  and  signal-to-noise  ratios. 
Specifically  we  presented  results  of  how  the  unknown  planet  radius  and  potentially 
poorly  known  observer-planet-star  phase  angle  can  impact  retrievals  of  atmospheric 
methane  abundance,  cloud  properties  and  surface  gravity,  among  others. 

Given  a  varying  planet  phase,  we  find  that  the  methane  abundance  can 
typically  only  be  constrained  to  the  correct  order  of  magnitude  at  SNR  of  20  or 
greater.  For  all  three  metallicities  considered  (lx,  lOx,  50x  solar),  only  the  case  with 
an  SNR  of  20  reproduces  to  the  correct  methane  abundance  at  all  phase  angles.  Low 
SNR  cases  merely  identify  the  presence  of  methane,  with  the  abundance  being  highly 
uncertain  across  several  orders  of  magnitude,  an  important  result  for  the  design  of 
future  space-based  missions  such  as  WFIRST. 
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Surface  gravity  appears  essentially  unconstrained.  The  top  cloud  in  a  two- 
layer  cloud  model  is  well  constrained,  but  is  indeterminate  on  the  pressure  gap 
between  cloud  layers,  indicating  that  a  one-cloud  model  might  be  better  suited  to  the 
examples  in  this  paper.  However  our  MCMC  methods  are  able  to  return  a  solution  for 
planet  radius  within  a  factor  of  two  of  the  true  value,  even  at  low  SNR  values. 
Surprisingly,  the  confidence  interval  on  the  radius  solution  decreases  with  increasing 
phase  angle.  Since  the  atmosphere  is  highly  forward  scattering  and  molecular  bands 
are  weak  at  more  crescent  phases,  clouds  become  less  important.  Retrievals  for  radius 
are  consequently  best  at  the  highest  phase  angles. 

We  find  that  knowledge  of  the  phase  angle,  and  therefore  its  elimination  as  a 
free  parameter,  does  not  significantly  improve  estimates  for  methane  abundance, 
cloud  parameters  or  gravity.  However  it  does  improve  the  radius  determination  by  a 
factor  of  two.  On  the  other  hand,  with  no  prior  knowledge  of  orbital  parameters,  the 
most  we  can  confidently  tell  about  the  phase  angle  from  retrieved  results  is  whether  it 
is  high  or  low.  Observations  with  comparable  SNR  and  multiple  phase  angles  can 
therefore  be  extremely  valuable  in  determining  both  orbital  characteristics,  if 
unknown  or  uncertain,  and  narrowing  down  planet  radius. 

Finally,  we  find  that  simultaneous  retrievals  on  observations  taken  at  multiple 
phase  angles  holds  promise  for  narrowing  the  solution  space  of  best-fit  models.  We 
estimate  this  using  an  intersection  criterion  and  find  a  steadily  increasing 
improvement  in  estimates  of  all  parameters,  even  for  generally  indeterminate 
parameters  such  as  surface  gravity  and  retrieved  phase  angle.  For  low  SNR  cases. 
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estimates  for  methane  abundance  can  be  improved  by  as  much  as  1-2  orders  of 
magnitude  if  multiple  observations  at  different  phase  angles  are  gathered,  a  fact  of 
interest  when  planning  future  space-based  observational  campaigns.  However,  it  is 
important  not  to  assign  too  much  importance  to  this,  since  even  though  bounds  on  the 
solution  may  decrease,  this  does  not  guarantee  the  accuracy  of  the  solution.  At  low 
SNRs,  the  recovered  methane  solution  is  far  separated  from  the  true  value  at  all  phase 
angles.  The  best  that  multiple  observations  at  low  SNR  can  hope  to  accomplish  is  the 
information  content  of  one  observation  at  high  SNR. 

Our  group  is  continuing  to  pursue  MCMC  methods  for  application  to 
reflected-light  spectral  data  in  the  context  of  a  wide  range  of  future  missions, 
including  WFIRST.  Future  work  will  focus  on  a  continued  improved  treatment  of 
clouds  and  hazes,  as  well  as  Raman  scattering,  although  we  expect  this  latter  effect  to 
be  minimal,  since  Raman  scattering  features  are  weak  in  the  visible  wavelengths 
[Karkoschka,  1994].  We  will  also  pursue  retrievals  to  determine  the  planetary 
temperature-pressure  profile  via  the  reflection  spectrum.  Finally,  we  will  perform 
simultaneous  retrievals  on  observations  taken  at  multiple  phase  angles,  an 
improvement  on  the  intersection  criterion  approximation  investigated  here. 
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5.7  Figures 
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Figure  5.1.  Illustration  of  the  degenerate  relationship  between  decreasing  planet  phase 
(increasing  phase  angle  a)  and  increasing  planet  radius,  in  yielding  an  equivalent  scattered  flux. 
If  the  planet  phase  is  unknown  or  uncertain,  a  larger  planet  at  a  crescent  phase  may  reflect 
essentially  the  same  amount  of  light  as  a  smaller  planet  at  a  fuller  phase. 
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Figure  5.2.  Illustrative  representation  of  two-layer  cloud  model  employed,  after  Marley  et  al. 
(2014). 
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Figure  5.3.  Candidate  target  planets  favorable  for  characterization  by  WFIRST.  The  M  sin  i  of 
each  planet,  as  determined  by  radial-velocity  measurements,  is  plotted  against  the  planet’s 
estimated  effective  temperature,  accounting  for  both  absorption  of  incident  flux  and  thermal 
evolution  as  described  in  Marley  et  al.  (2014).  Color  handing  indicates  approximate  effective 
temperature  regimes  where  various  clouds  will  dominate  the  reflected  flux  signal.  Dashed  boxes 
highlight  planets  discussed  in  this  work.  Jupiter  is  indicated  by  the  gold  dot.  The  green  dot 
indicates  the  effective  temperature  (but  not  the  mass,  which  is  lower)  of  Uranus. 


Wavelength  (microns) 

Figure  5.4.  Model  noise-free  contrast  spectra  for  the  two  test  cases.  Planet  A  (red)  and  HD 
192310c  (50x  metallicity,  blue),  for  a  spectral  resolution  of  R  =  70.  Dashed  spectra  indicate  an 
observation  at  quadrature;  methane  absorption  features  are  still  notable. 
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Figure  5.5.  Contrast  spectra  for  Planet  A  at  SNR  values  of  (a)  5  and  (b)  20.  Spectra  are 
generated  at  a  spectral  resolution  of  R  =  70  and  a  phase  angle  of  30°.  Red  represents  the  truth 
spectra  and  blue  error  bars  represent  notional  instrument  noise  during  observation  [Robinson  et 
a/.,  2015]. 
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Figure  5.6.  Contrast  spectra  for  HD  192310c  generated  at  metallicity  values  of  (a)  lx,  (b)  lOx  and 
(c)  50x  that  of  the  Sun.  The  difference  in  the  methane  absorption  signature  at  -0.9  micron  is 
evident.  Spectra  are  generated  at  phase  angles  of  0°,  SNR  =  20  and  spectral  resolution  of  R  =  70. 
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Figure  5.7.  Best-fit  contrast  spectra  for  Planet  A  at  SNR  =  10  and  varying  phase  angles:  (a)  30°; 
(b)  60°;  (c)  90°;  (d)  120°.  Model  spectra  are  calculated  from  the  best  198,000  samples  (1500 
iterations)  of  the  330,000  final  sample  chain  (2500  iterations).  The  median  spectrum  (blue) 
matches  well  to  the  truth  spectrum  (black).  16-84%  (dark  red)  and  4.5-95.5%  (light  red) 
percentile  range  of  recovered  solutions  are  also  shown.  Good  matches  to  the  model  truth  spectra 
are  seen  in  all  cases,  even  at  relatively  low  contrast  signals  for  large  phase  angles.  Note  the 
differing  (smaller)  vertical  scales  between  (a)  and  (d). 
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Contrast  Contrast  Contrast 


Figure  5.8.  Best-fit  contrast  spectra  for  HD  192310c  at  a  metallicity  of  lx  solar  and  varying  SNR 
and  phase  angle:  (a)  SNR  =  5,  a  =  30°;  (b)  SNR  =  5,  a  =  90°;  (c)  SNR  =  10,  a  =  30°;  (d)  SNR  =  10, 
a  =  90°;  (e)  SNR  =  20,  a  =  30°;  (f)  SNR  =  20,  a  =  90°.  Note  the  lack  of  the  characteristic  methane 
absorption  signal  at  0.9  microns  for  high  phase  angles.  Description  of  colors  is  as  in  Figure  5.7. 
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Figure  5.9.  Best-fit  contrast  spectra  for  HD  192310c  at  a  metallicity  of  lOx  solar  and  varying 
SNR  and  phase  angle:  (a)  SNR  =  5,  a  =  30°;  (b)  SNR  =  5,  a  =  90°;  (c)  SNR  =  10,  a  =  30°;  (d)  SNR 
=  10,  a  =  90°;  (e)  SNR  =  20,  a  =  30°;  (f)  SNR  =  20,  a  =  90°.  Good  matches  to  the  model  truth 
spectra  are  seen  in  all  cases,  even  at  relatively  low  contrast  signals  for  large  phase  angles.  The 
improvement  in  recovered  signal  is  evident  with  increasing  SNR  from  (a)  through  (f). 
Description  of  colors  is  as  in  Figure  5.7. 
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Figure  5.10.  Best-fit  contrast  spectra  for  HD  192310c  at  a  metallicity  of  50x  solar  and  varying 
SNR  and  phase  angle:  (a)  SNR  =  5,  a  =  30°;  (b)  SNR  =  5,  a  =  90°;  (c)  SNR  =  10,  a  =  30°;  (d)  SNR 
=  10,  a  =  90°;  (e)  SNR  =  20,  a  =  30°;  (f)  SNR  =  20,  a  =  90°.  Good  matches  to  the  model  truth 
spectra  are  seen  in  all  cases,  even  at  relatively  low  contrast  signals  for  large  phase  angles.  The 
improvement  in  recovered  signal  is  evident  with  increasing  SNR  from  (a)  through  (f). 
Description  of  colors  is  as  in  Figure  5.7. 
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Figure  5.11.  Posterior  probability  distribution  plot  for  all  eleven  parameters  retrieved  by  the 
MCMC  algorithm  (Table  5.1),  for  the  case  of  Planet  A  at  SNR  =  10  and  phase  angle  60°  (Figure 
5.7c).  Darker  regions  represent  higher  probability.  Blue  dots  represent  true  values.  The 
distributions  are  drawn  from  all  remaining  samples  after  the  MCMC  burn-in  chains  (first  1000 
chains)  are  discarded.  The  diagonal  of  the  plot  represents  the  marginalized  probability 
distributions  for  each  parameter.  Best-fit  values  in  log  space  (except  for  phase  angle)  are  shown. 
See  text  for  references  to  text  markers  a-c.  The  error  bars  indicate  (left  to  right)  the  16%,  84% 
and  50%  quantiles,  i.e.,  this  is  the  68%  confidence  interval. 
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Figure  5.12.  Highlight  of  pertinent  parameter  relationships  with  planet  phase  angle,  excerpted 
from  Figure  5.11.  The  top  row  shows  probability  distributions  of  phase  angle  against  methane 
abundance,  surface  gravity  and  planet  radius.  The  bottom  row  shows  distributions  of  phase 
angle  against,  respectively,  the  pressure  of  the  bottom  cloud  (P2)  in  the  two-layer  model  by 
Marley  et  al.  (2014),  the  pressure  difference  between  the  two  cloud  layers  (dPi)  and  the  pressure 
difference  across  the  top  cloud  layer  (dP2).  All  parameters  are  in  log  space  (except  phase  angle). 
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Figure  5.13.  The  relationship  of  planet  radius  with  changing  planet  phase  angle,  for  the  Planet  A 
test  case,  for  an  SNR  of  20  and  a  truth  phase  angle  (blue  square)  of:  (a)  30°;  (b)  60°,  (c)  90°  and 
(d)  120°.  The  MCMC  algorithm  retrieves  phase  angle  and  radius  values  close  to  the  true  value 
for  planet  phases  close  to  full  phase  and  past  quadrature.  Best-fit  values  for  each  case  are 
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indicated  above  the  figure;  superscripts  and  subscripts  to  this  value  indicate  upper  and  lower 
bounds  returned  from  the  posterior  probability  plot  diagonals  (e.g.  Figure  5.11). 
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Figure  5.14.  Summary  of  all  results  for  Planet  A  test  case.  True  phase  angle  varies  from  0°  to 
120°.  SNR  varies  between  5  (red),  10  (green)  and  20  (blue).  Error  bars  enclose  the  68% 
confidence  interval  as  defined  by  MCMC  posterior  probability  distributions  similar  to  Figure 
5.11.  Solid  dots  denote  the  best-fit  values.  A  black  dashed  line  denotes  true  values  from  the 
Planet  A  model. 
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Figure  5.15.  Summary  of  all  results  for  the  HD  192310c  test  case  with  lx  metallicity  of  the  Sun. 
True  phase  angle  varies  from  0°  to  120°.  SNR  varies  between  5  (red),  10  (green)  and  20  (blue). 
Error  bars  are  as  in  Figure  5.14.  Solid  dots  denote  the  best-fit  values.  A  black  dashed  line 
denotes  true  values  from  the  HD  192310c  model. 
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Figure  5.16.  Summary  of  all  results  for  the  HD  192310c  test  case  with  lOx  metallicity  of  the  Sun. 
True  phase  angle  varies  from  0°  to  120°.  SNR  varies  between  5  (red),  10  (green)  and  20  (blue). 
Error  bars  and  solid  dots  are  as  in  Figure  5.14.  A  black  dashed  line  denotes  true  values  from  the 
HD  192310c  model. 


Figure  5.17.  Summary  of  all  results  for  the  HD  192310c  test  case  with  50x  metallicity  of  the  Sun. 
True  phase  angle  varies  from  0°  to  120°.  SNR  varies  between  5  (red),  10  (green)  and  20  (blue). 
Error  bars  and  solid  dots  are  as  in  Figure  5.14.  A  black  dashed  line  denotes  true  values  from  the 
HD  192310c  model. 
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Figure  5.18.  Results  for  the  lOx  metallicity  case  of  HD  192310c,  at  SNR  =  20,  for  an  unbounded 
phase  angle  case  (red)  and  a  bounded  case  (green),  where  for  the  bounded  case,  the  phase  angle 
can  vary  by  no  more  than  10°  from  the  true  value.  No  impact  to  cloud  property  or  methane 
abundance  retrievals  are  seen,  however,  the  retrieved  radius  of  the  planet  improves  by  a  factor 
of  two. 


Figure  5.19.  Illustration  of  the  intersection  criterion  for  surface  gravity  and  planet  radius  (note 
log  units  for  both).  Data  is  from  retrievals  for  HD  192310c,  metallicity  lOx  solar  and  an  SNR  of 
5.  Phase  angles  are  shown  in  the  legend.  The  dashed  black  line  indicates  true  values.  The 
highlighted  yellow  area  represents  the  region  of  solutions  common  to  both  retrievals  (intersection 
criterion);  it  can  be  seen  to  visibly  improve  error  bars  on  both  cumulative  solutions. 
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<—  Phase  angle  of  first  observation  (deg)  — > 

Figure  5.20.  Improvement  in  68%  confidence  interval  ranges  for  methane  abundance,  surface 
gravity,  planet  radius,  retrieved  phase  angle  and  two-cloud  top  pressures;  peaks  represent 
maximum  improvement  between  two  observations  taken  at  different  phase  angles.  Results  are 
for  HD  192310c  (lOx  metallicity  case)  at  SNR  =  5.  The  first  phase  angle  is  represented  on  the  x- 
axis;  the  second  as  colored  lines  (see  figure  legend).  Units  of  the  relative  improvement  are 
indicated  in  the  figure  titles. 
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Figure  5.21.  Improvement  in  68%  confidence  intervals  for  an  intersection  criterion  applied  to 
three  randomly  chosen  phase-varying  observations  of  (A)  Planet  A  and  (B)  HD  192310c,  lOx 
metallicity  case,  compared  to  a  randomly  chosen  single  observation.  SNR  values  of  5-20  are 
shown  (blue:  SNR  5,  green:  SNR  10,  red:  SNR  20).  See  text  for  references  to  text  markers  d-h. 
This  plot  may  be  used  to  determine  trends  for  improvement  in  uncertainty  estimates  by 
retrieving  against  multiple  datasets  collected  at  differing  phase  angles.  For  example,  for  gravity, 
most  cases  show  no  improvement  in  gravity  estimates  (improvement  clusters  around  0  ms'^), 
regardless  of  SNR,  when  confidence  intervals  from  multiple  observations  are  stacked  together, 
regardless  of  the  phase  angle  of  the  single  observation.  However,  for  methane  abundance,  a 
significant  number  of  SNR  5  cases  (blue)  show  at  0.5-1  order  of  magnitude  improvement  when 
the  intersection  criterion  is  applied  to  multiple  observations.  A  similar  trend  is  noted  for  higher 
SNR  cases,  although  the  number  of  cases  that  note  this  improvement  drops  off  as  expected. 
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Figure  5.22.  SNR  versus  phase  angle  relationships  for  all  four  test  cases,  with  respect  to  methane 
abundance.  Colors  indicate  the  size  of  68%  confidence  interval  error  bars  in  orders  of 
magnitude.  Values  originate  from  posterior  probability  distributions  similar  to  Figure  5.11.  Blue 
colors  indicate  a  tighter  confidence  interval,  i.e.,  a  better  retrieval. 
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Figure  5.23.  SNR  versus  phase  angle  relationships  for  all  four  test  cases,  with  respect  to  surface 
gravity.  Colors  indicate  the  size  of  68%  confidence  interval  error  bars  in  units  of  ms'^.  Values 
originate  from  posterior  probability  distributions  similar  to  Figure  5.11.  Blue  colors  indicate  a 
tighter  confidence  interval,  i.e.,  a  better  retrieval. 
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Figure  5.24.  SNR  versus  phase  angle  relationships  for  all  four  test  cases,  with  respect  to  planet 
radius.  Colors  indicate  the  size  of  68%  confidence  interval  error  bars  in  units  of  Rj.  Values 
originate  from  posterior  probability  distributions  similar  to  Figure  5.11.  Blue  colors  indicate  a 
tighter  confidence  interval,  i.e.,  a  better  retrieval. 
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Figure  5.25.  SNR  versus  phase  angle  relationships  for  all  four  test  cases,  with  respect  to  inferred 
phase  angle.  Colors  indicate  the  size  of  68%  confidence  interval  error  bars  in  units  of  degrees. 
Values  originate  from  posterior  probability  distributions  similar  to  Figure  5.11.  Blue  colors 
indicate  a  tighter  confidence  interval,  i.e.,  a  better  retrieval. 
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Figure  5.26.  SNR  versus  phase  angle  relationships  for  all  four  test  cases,  with  respect  to  the 
pressure  at  the  top  of  the  bottom  cloud.  Colors  show  size  of  68%  confidence  interval  error  bars 
in  orders  of  magnitude.  Values  originate  from  posterior  probability  distributions  similar  to 
Figure  5.11.  Blue  colors  indicate  a  tighter  confidence  interval,  i.e.,  a  better  retrieval. 
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Figure  5.27.  SNR  versus  phase  angle  relationships  for  all  four  test  cases,  with  respect  to  the 
pressure  at  the  top  of  the  top  cloud.  Colors  show  size  of  68%  confidence  interval  error  bars  in 
orders  of  magnitude.  Values  originate  from  posterior  probability  distributions  similar  to  Figure 
5.11.  Blue  colors  indicate  a  tighter  confidence  interval,  i.e.,  a  better  retrieval. 
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